Blue Genes Analysis

# Read in libraries 
library(here); library(tidyverse); library(emmeans);
library(lme4); library(lmerTest); library(kableExtra);
library(car); library(sjPlot); library(MuMIn); library(EnvStats);
library(merTools); library(patchwork); library(RColorBrewer)

# Read in compiled trait data
source(here("supp_code", "CompileTraitData.R"))

Data Subsetting

For most analyses we are going to need to subset the data into the competition data set (all Corn Island – these were the only genotypes that had competition; and the no competition pots). For the no competition pots, we are also dropping out the Blackwater genotypes because they can’t really be separated into age and location cohorts. We are also going to focus in on levels 1-4 where most clones survived.

Calculate the overall survival rate by flooding level

level	n	extinct	survival
1	112	7	93.8%
2	112	6	94.6%
3	112	6	94.6%
4	112	22	80.4%
5	110	80	27.3%
6	112	109	2.7%

Summarizing by cohort and genotype for the all Corn Island dataset

age	location	n
ancestral	corn	160
modern	corn	162

genotype	n
ca9	6
ca8	12
ca2	13
cm8	13
ca3	14
cm6	14
cm7	14
ca6	15
cm3	15
cm5	15
ca1	16
ca5	16
ca7	16
cm1	16
cm2	16
cm9	16
cm10	20
cm4	23
ca10	26
ca4	26

Summarizing by cohort and genotype for the no competition dataset

age	location	n
ancestral	corn	79
ancestral	kirkpatrick	44
modern	corn	84
modern	kirkpatrick	48

genotype	n
ca9	1
ca8	3
ca3	4
cm8	5
cm9	5
km4	5
ka1	6
ka2	6
ka3	6
ca2	7
km2	7
ca7	8
cm2	8
cm3	8
cm6	8
ka4	8
km5	8
km6	8
ca1	9
ca6	9
cm1	9
ka5	9
ka6	9
km1	9
ca5	10
cm4	10
cm5	10
cm7	10
cm10	11
km3	11
ca4	12
ca10	16

Trait analysis for no competition pots

For the overall trait analysis, we are going to do two separate analyses, one for each data set. Because there are so many potential interactions between factors (environmental or genetic), we will use a backward selection model approach to see what interactions are important in explaining variation in the data. Here is the general approach:

• Fit the most complex model with respect to fixed effects
  • 5-way interaction between \(\text{CO}_2\) (2 levels: ambient, elevated), salinity (2 levels: GCREW, freshwater), elevation (continuous), age cohort (2 levels: ancestral, descendant), and provenance (2 levels: Kirkpatrick Marsh, Corn Island)
  • quadratic term for elevation (no interactions)
  • initial weight of the propagule, temporal planting group (2 levels: late June, early July), lab of origin (2 levels: ND, UTK)
• In this model also include random intercepts for genotype and for chamber
  • For some models, the variance explained by chamber is basically zero and can be dropped from the model
• Use lmerTest::step for perform backwards stepwise selection on the fixed effects of that model
  • Use the argument keep to force the stepwise selection procedure to keep the additive fixed effects of \(\text{CO}_2\), salinity, elevation, age cohort, and provenance because these are things that we purposefully manipulated
  • This means that the only terms that can be dropped are interactions and covariates related to the experimental set-up such as initial weight of the propagule
• Given the final model, then test to see if there are genotype by environment interactions that can explain significantly more variation
  • Do this by fitting models with random slopes using terms (up to two-way interactions) that were kept in the model for \(\text{CO}_2\), salinity, and elevation
  • Test to see if the genotype by environment should be kept using a significance test with anova() while setting the argument refit = FALSE to fit models using restricted maximum likelihood for these comparsions

Aboveground biomass

# Center and scale continuous variables and subset to just aboveground biomass
# greater than 0
bg_nocomp %>% 
  mutate(elevation_sc = scale(elevation)[,1],
         weight_init_sc = scale(weight_init)[,1]) %>% 
  filter(extinct == 0)-> bg_nocomp

# Fit the most complex model -- sqrt of agb biomass seems like the best
# transformation
agb_mod <- lmer(sqrt(agb_scam) ~ weight_init_sc + date_planted_grp + origin_lab +
                  (age + location + co2 + salinity + elevation_sc)^5 + I(elevation_sc^2) +
                  (1|site_frame) + (1|genotype), data = bg_nocomp)

## boundary (singular) fit: see ?isSingular

# Warning indicates that one of the random intercepts is estimated to be zero
VarCorr(agb_mod) # It's chamber, so we'll drop it

##  Groups     Name        Std.Dev.
##  genotype   (Intercept) 0.24224 
##  site_frame (Intercept) 0.00000 
##  Residual               0.51650

# Fit a model without that term
agb_mod2 <- update(agb_mod, .~.-(1|site_frame))

# Do a stepwise model selection
# step(agb_mod2, keep = c("co2", "elevation_sc", "salinity", "age", "location"))

# There's a model fitting issue because the model is too complicated. Let's see
# if there's a significant 5-way interaction
Anova(agb_mod2)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: sqrt(agb_scam)
##                                          Chisq Df Pr(>Chisq)    
## weight_init_sc                          4.5169  1   0.033562 *  
## date_planted_grp                        1.1259  1   0.288661    
## origin_lab                              0.0399  1   0.841596    
## age                                     0.4286  1   0.512701    
## location                                0.2874  1   0.591898    
## co2                                     4.2801  1   0.038561 *  
## salinity                                0.3342  1   0.563194    
## elevation_sc                           68.4369  1  < 2.2e-16 ***
## I(elevation_sc^2)                      25.5516  1  4.307e-07 ***
## age:location                            1.5009  1   0.220539    
## age:co2                                 0.0240  1   0.876871    
## age:salinity                            0.1001  1   0.751663    
## age:elevation_sc                        0.0282  1   0.866728    
## location:co2                            0.6727  1   0.412121    
## location:salinity                       0.0117  1   0.913832    
## location:elevation_sc                   0.5459  1   0.459995    
## co2:salinity                            0.7830  1   0.376220    
## co2:elevation_sc                        4.1267  1   0.042212 *  
## salinity:elevation_sc                   0.0245  1   0.875622    
## age:location:co2                        2.4417  1   0.118148    
## age:location:salinity                   0.0659  1   0.797354    
## age:location:elevation_sc               0.0807  1   0.776346    
## age:co2:salinity                        4.9315  1   0.026371 *  
## age:co2:elevation_sc                    0.0013  1   0.971258    
## age:salinity:elevation_sc               0.0832  1   0.772974    
## location:co2:salinity                   0.5042  1   0.477659    
## location:co2:elevation_sc               0.1599  1   0.689287    
## location:salinity:elevation_sc          1.3568  1   0.244085    
## co2:salinity:elevation_sc               0.0737  1   0.785970    
## age:location:co2:salinity               0.3398  1   0.559959    
## age:location:co2:elevation_sc           0.2862  1   0.592665    
## age:location:salinity:elevation_sc      0.4013  1   0.526431    
## age:co2:salinity:elevation_sc           7.6946  1   0.005539 ** 
## location:co2:salinity:elevation_sc      2.2696  1   0.131935    
## age:location:co2:salinity:elevation_sc  0.0104  1   0.918683    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Nope! Drop it and try stepwise again
agb_mod3 <- lmer(sqrt(agb_scam) ~ weight_init_sc + date_planted_grp + origin_lab +
                  (age + location + co2 + salinity + elevation_sc)^4 + I(elevation_sc^2) +
                  (1|genotype), data = bg_nocomp)

step(agb_mod3, keep = c("co2", "elevation_sc", "salinity", "age", "location"))

## Backward reduced random-effect table:
## 
##                Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)   
## <none>                      37 -207.59 489.19                        
## (1 | genotype)          0   36 -212.58 497.17 9.9812  1   0.001582 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                                    Eliminated Sum Sq Mean Sq NumDF   DenDF
## origin_lab                                  1 0.0103  0.0103     1  31.802
## age:location:co2:elevation_sc               2 0.0738  0.0738     1 183.667
## age:location:co2:salinity                   3 0.0958  0.0958     1 177.380
## age:location:salinity:elevation_sc          4 0.1162  0.1162     1 181.096
## age:location:salinity                       5 0.0155  0.0155     1 184.257
## age:location:elevation_sc                   6 0.0170  0.0170     1 182.359
## date_planted_grp                            7 0.2300  0.2300     1 190.751
## location:co2:salinity:elevation_sc          8 0.5739  0.5739     1 201.997
## location:co2:elevation_sc                   9 0.0440  0.0440     1 192.385
## location:co2:salinity                      10 0.2111  0.2111     1 184.830
## location:salinity:elevation_sc             11 0.3550  0.3550     1 189.972
## location:salinity                          12 0.0194  0.0194     1 192.593
## location:elevation_sc                      13 0.1527  0.1527     1 189.981
## age:location:co2                           14 0.6496  0.6496     1 196.121
## location:co2                               15 0.2865  0.2865     1 196.906
## age:location                               16 0.5040  0.5040     1  24.407
## weight_init_sc                              0 1.7330  1.7330     1 207.800
## location                                    0 0.0558  0.0558     1  25.858
## I(elevation_sc^2)                           0 6.4346  6.4346     1 198.512
## age:co2:salinity:elevation_sc               0 1.8198  1.8198     1 206.044
##                                    F value    Pr(>F)    
## origin_lab                          0.0386  0.845437    
## age:location:co2:elevation_sc       0.2781  0.598571    
## age:location:co2:salinity           0.3622  0.548069    
## age:location:salinity:elevation_sc  0.4412  0.507405    
## age:location:salinity               0.0588  0.808599    
## age:location:elevation_sc           0.0651  0.798966    
## date_planted_grp                    0.8840  0.348286    
## location:co2:salinity:elevation_sc  2.2112  0.138573    
## location:co2:elevation_sc           0.1666  0.683634    
## location:co2:salinity               0.8054  0.370644    
## location:salinity:elevation_sc      1.3544  0.245970    
## location:salinity                   0.0740  0.785875    
## location:elevation_sc               0.5838  0.445767    
## age:location:co2                    2.4948  0.115837    
## location:co2                        1.0872  0.298363    
## age:location                        1.9062  0.179896    
## weight_init_sc                      6.5577  0.011153 *  
## location                            0.2110  0.649789    
## I(elevation_sc^2)                  24.3478 1.698e-06 ***
## age:co2:salinity:elevation_sc       6.8859  0.009338 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## sqrt(agb_scam) ~ weight_init_sc + age + location + co2 + salinity + 
##     elevation_sc + I(elevation_sc^2) + (1 | genotype) + age:co2 + 
##     age:salinity + age:elevation_sc + co2:salinity + co2:elevation_sc + 
##     salinity:elevation_sc + age:co2:salinity + age:co2:elevation_sc + 
##     age:salinity:elevation_sc + co2:salinity:elevation_sc + age:co2:salinity:elevation_sc

# Fit the proposed stepwise model
agb_mod4 <- lmer(sqrt(agb_scam) ~ weight_init + I(elevation^2) + location +
    age*co2*salinity*elevation + (1 | genotype), data = bg_nocomp)

# See if there are random effects models that are an improvement of this model
agb_mod5 <- update(agb_mod4, .~.-(1|genotype) + (co2*elevation_sc|genotype))

## boundary (singular) fit: see ?isSingular

agb_mod6 <- update(agb_mod4, .~.-(1|genotype) + (co2*salinity|genotype))

## boundary (singular) fit: see ?isSingular

agb_mod7 <- update(agb_mod4, .~.-(1|genotype) + (salinity*elevation_sc|genotype))

## boundary (singular) fit: see ?isSingular

agb_mod8 <- update(agb_mod4, .~.-(1|genotype) + (co2+elevation_sc|genotype))

## boundary (singular) fit: see ?isSingular

agb_mod9 <- update(agb_mod4, .~.-(1|genotype) + (co2+salinity|genotype))

## boundary (singular) fit: see ?isSingular

agb_mod10 <- update(agb_mod4, .~.-(1|genotype) + (salinity+elevation_sc|genotype))

## boundary (singular) fit: see ?isSingular

agb_mod11 <- update(agb_mod4, .~.-(1|genotype) + (co2|genotype)) 

## boundary (singular) fit: see ?isSingular

agb_mod12 <- update(agb_mod4, .~.-(1|genotype) + (salinity|genotype))

## boundary (singular) fit: see ?isSingular

agb_mod13 <- update(agb_mod4, .~.-(1|genotype) + (elevation_sc|genotype))

## boundary (singular) fit: see ?isSingular

# So final model is the same as the stepwise proposed model
agb_final <- agb_mod4

# Look at model diagnostics
plot_model(agb_final, type = "diag")

## [[1]]

## `geom_smooth()` using formula 'y ~ x'

## 
## [[2]]
## [[2]]$genotype

## `geom_smooth()` using formula 'y ~ x'

## 
## 
## [[3]]

## 
## [[4]]

## `geom_smooth()` using formula 'y ~ x'

# Look at random effects
plot_model(agb_final, type = "re")

# Refit model and get Type 3 Anova table
Anova(agb_final, type = 3)

## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: sqrt(agb_scam)
##                              Chisq Df Pr(>Chisq)    
## (Intercept)                 9.1719  1  0.0024576 ** 
## weight_init                 6.5577  1  0.0104433 *  
## I(elevation^2)             24.3478  1  8.042e-07 ***
## location                    0.2110  1  0.6459501    
## age                         2.0667  1  0.1505426    
## co2                         0.0603  1  0.8059773    
## salinity                    2.4430  1  0.1180472    
## elevation                  11.5307  1  0.0006845 ***
## age:co2                     5.8183  1  0.0158601 *  
## age:salinity                3.7727  1  0.0520955 .  
## co2:salinity                3.5302  1  0.0602594 .  
## age:elevation               1.8766  1  0.1707246    
## co2:elevation               0.1410  1  0.7072758    
## salinity:elevation          0.9067  1  0.3409938    
## age:co2:salinity            9.8915  1  0.0016604 ** 
## age:co2:elevation           4.3381  1  0.0372683 *  
## age:salinity:elevation      2.2560  1  0.1330943    
## co2:salinity:elevation      2.9289  1  0.0870051 .  
## age:co2:salinity:elevation  6.8859  1  0.0086878 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Extract predicted means for model
agb_means <- emmeans::emmeans(agb_final, ~co2:salinity:age:elevation, at= list(elevation = seq(0.15, 0.55, 0.01)), type = "response")

agb_means_out <- summary(agb_means)
agb_means_out %>% 
  mutate(salinity = case_when(salinity == "fresh" ~ "Freshwater Site (6ppt)",
                              T ~ "Brackish Site (8ppt)"),
         age = case_when(age == "modern" ~ "Descendant cohort (2000-2020)",
                         T ~ "Ancestral cohort (1900-1960)")) -> agb_means_out

bg_nocomp_plot <- bg_nocomp %>% 
  mutate(salinity = case_when(salinity == "fresh" ~ "Freshwater Site (6ppt)",
                              T ~ "Brackish Site (8ppt)"),
         age = case_when(age == "modern" ~ "Descendant cohort (2000-2020)",
                         T ~ "Ancestral cohort (1900-1960)"))


agb_means_out %>% 
  ggplot(aes(x = elevation, y = response, color = co2)) + 
  geom_line(size = 1.5) +
  #geom_line(aes(x = elevation, y = lower.CL), linetype = "dashed") + 
  #geom_line(aes(x = elevation, y = upper.CL), linetype = "dashed") + 
  facet_grid(age~salinity) +
  geom_point(data = bg_nocomp_plot, aes(x = elevation, y = agb_scam), alpha = 0.5) +
  theme_bw() +
  scale_color_manual(values = c("gray47", "#66a61e")) +
  ylab("Aboveground biomass (g)") +
  xlab("Elevation (m NAVD88)") +
  guides(color=guide_legend(title=expression(CO[2]))) +
  theme(legend.position = "top")-> a

# Comparison to a null model that does not include age, location or genotype
agb_null <- lm(sqrt(agb_scam) ~ weight_init_sc +
                  (co2 + salinity + elevation_sc)^3 + I(elevation_sc^2), data = bg_nocomp)

# Get R2 from model
agb_null_r2 <- r.squaredGLMM(agb_null)[2]

## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.

# Compare to conditional R2 from full model
agb_r2 <- r.squaredGLMM(agb_final)[2]

# Create enough colors for each genotype to have one
n <- length(unique(bg_nocomp$genotype))
qual_col_pals = brewer.pal.info[brewer.pal.info$category == 'qual',]
col_vector = unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals)))

# Create caterpillar graph
REsim(agb_final) %>%
  mutate(groupID = fct_reorder(groupID, mean)) %>% 
  ggplot(aes(x = groupID)) +
  geom_pointrange(aes(ymin = mean - 2*sd, ymax = mean + 2*sd, xmin = groupID, xmax = groupID, y = mean, color = groupID)) +
  coord_flip() +
  ylab("Random effect") +
  xlab("Genotype") +
  geom_hline(aes(yintercept = 0), linetype = "dashed") +
  theme_bw() + scale_x_discrete(position = "top")  +
  scale_color_manual(values = col_vector) +
  theme(legend.position = "none")-> b
  #theme(axis.text.y=element_blank(),
        #axis.ticks.y = element_blank()) -> b


cowplot::plot_grid(a,b, rel_widths = c(4,1)) -> agb_plot

Root:shoot ratio

# Calculate root-to-shoot ratio on no comp pots
bg_nocomp %>% 
  mutate(rs = total_bg / agb_scam) -> bg_nocomp

# Plot and see if there are outliers
bg_nocomp %>% 
  ggplot(aes(x = rs)) +
  geom_histogram() +
  geom_rug() +
  xlab("root-to-shoot ratio")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 2 rows containing non-finite values (stat_bin).

# There are two pots that have really high r:s and one that seems a bit high.
# Take a look at these further.
bg_nocomp %>% 
  filter(rs > 4.5)

##   pot_no      beta bg_scam10 bg_sppa10 total_bg site_frame level position
## 1     98 0.9392815     1.460        NA    5.552    fresh_3     1        2
## 2    314        NA     0.084        NA    0.084    fresh_7     4        2
## 3   1936 0.9584345     0.151        NA    0.371    gcrew_6     4        6
##       co2 salinity comp elevation genotype location depth_seed         cohort
## 1 ambient    fresh    0 0.5260000      cm8     corn          0    corn_modern
## 2 ambient    fresh    0 0.1980000      cm3     corn          0    corn_modern
## 3 ambient     salt    0 0.1853333      ca1     corn      20.75 corn_ancestral
##   origin_lab date_cloned_grp date_planted_grp weight_init width_init nodes_init
## 1         nd               1                1       0.667       2.29          2
## 2        utk               1                1       0.724       0.70          1
## 3        utk               1                2       0.390       1.44          1
##   agb_scam agb_sppa dens_scam_live dens_scam_dead height_sppa height_scam_tot
## 1    1.186        0             12              0          NA        40.08333
## 2    0.007        0              1              0          NA              NA
## 3    0.047        0              1              0          NA        35.50000
##   height_scam_grn width_scam_mid width_scam_bas prop_twist soil_to_pot
## 1        38.16667       1.583333          1.775  0.4166667         6.0
## 2              NA             NA             NA         NA         2.5
## 3        34.00000       1.900000          2.300  0.0000000         5.0
##            id extinct       age elevation_sc weight_init_sc        rs
## 1 fresh_3_1_2       0    modern     1.392848     -0.7175561  4.681282
## 2 fresh_7_4_2       0    modern    -1.274990     -0.6410989 12.000000
## 3 gcrew_6_4_6       0 ancestral    -1.378016     -1.0891112  7.893617

# It seems like pot_no 98 is real -- just a lot of bg compared to ag (on level
# 1)

# The other two seem to be more of experimental artifacts (i.e. initial
# propagule weight is driving total belowground weight at the end), so we drop
# those two

bg_nocomp %>% 
  filter(rs < 6) -> bg_nocomp_rs

# There also appear to be some that are NAs
bg_nocomp %>% 
  filter(!complete.cases(rs))

##   pot_no beta bg_scam10 bg_sppa10 total_bg site_frame level position      co2
## 1    205   NA        NA        NA       NA    fresh_5     2        5 elevated
## 2   1700   NA        NA        NA       NA    gcrew_4     3        4 elevated
##   salinity comp elevation genotype location depth_seed         cohort
## 1    fresh    0 0.3973333      ca2     corn      18.75 corn_ancestral
## 2     salt    0 0.2871667      ca4     corn      20.75 corn_ancestral
##   origin_lab date_cloned_grp date_planted_grp weight_init width_init nodes_init
## 1        utk               1                1       0.733       1.71          1
## 2        utk               1                2       0.656       1.22          2
##   agb_scam agb_sppa dens_scam_live dens_scam_dead height_sppa height_scam_tot
## 1 7.287000        0             32              2          NA        47.51562
## 2 8.005503        0             27              3          NA        48.67308
##   height_scam_grn width_scam_mid width_scam_bas prop_twist soil_to_pot
## 1        45.42188       3.096875       3.287500  0.4375000         4.5
## 2        47.13462       3.326923       3.896154  0.3461538         4.0
##            id extinct       age elevation_sc weight_init_sc rs
## 1 fresh_5_2_5       0 ancestral    0.3463181     -0.6290267 NA
## 2 gcrew_4_3_4       0 ancestral   -0.5497391     -0.7323110 NA

# Need to figure this out on data checking!!! I think these were pots with
# messed up bg processing. Move forward with this for now, but I'll go back and
# fix these later

# Change names of salinity factor for presentation
bg_nocomp_rs %>% 
  mutate(salinity = case_when(salinity == "fresh" ~ "Freshwater Site (6ppt)",
                              T ~ "Brackish Site (8ppt)")) -> bg_nocomp_rs

# Scale elevation
bg_nocomp_rs$elevation_sc <- scale(bg_nocomp_rs$elevation)[,1]
centering <- attr(scale(bg_nocomp_rs$elevation), "scaled:center")
scaling <- attr(scale(bg_nocomp_rs$elevation), "scaled:scale")

## Fit model and model selection (fixed effects) ####

# Most complex model -- log is the best transformation
rs_mod <- lmer(log(rs) ~ weight_init + date_planted_grp + origin_lab +
                  (age + location + co2 + salinity + elevation)^5 + I(elevation^2) +
                  + (1|site_frame) + (1|genotype), data = bg_nocomp_rs)

# Do stepwise model selection
# step(rs_mod, keep = c("co2", "salinity", "elevation", "age", "location"))

# Error because model is too complex -- Try with just 4-way interactions
# Most complex model -- log is the best transformation
rs_mod2 <- lmer(log(rs) ~ weight_init + date_planted_grp + origin_lab +
                  (age + location + co2 + salinity + elevation)^4 + I(elevation^2) +
                  + (1|site_frame) + (1|genotype), data = bg_nocomp_rs)

# Do stepwise model selection
# step(rs_mod2, keep = c("co2", "salinity", "elevation", "age", "location")) # Still mad

Anova(rs_mod2) # Could drop down to 3-way interactions

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: log(rs)
##                                   Chisq Df Pr(>Chisq)    
## weight_init                      0.4123  1  0.5207998    
## date_planted_grp                 0.9076  1  0.3407395    
## origin_lab                       0.2156  1  0.6424063    
## age                              0.0785  1  0.7792791    
## location                         0.9893  1  0.3199205    
## co2                              6.5209  1  0.0106614 *  
## salinity                         1.5615  1  0.2114511    
## elevation                       17.5362  1  2.819e-05 ***
## I(elevation^2)                   3.5973  1  0.0578750 .  
## age:location                     0.1646  1  0.6849751    
## age:co2                          0.0388  1  0.8438226    
## age:salinity                     0.5912  1  0.4419611    
## age:elevation                    1.7376  1  0.1874398    
## location:co2                     0.0051  1  0.9432056    
## location:salinity                2.9562  1  0.0855463 .  
## location:elevation               5.1197  1  0.0236561 *  
## co2:salinity                     0.5831  1  0.4451132    
## co2:elevation                    0.0537  1  0.8166909    
## salinity:elevation              13.8903  1  0.0001938 ***
## age:location:co2                 0.6379  1  0.4244705    
## age:location:salinity            1.6866  1  0.1940492    
## age:location:elevation           0.5785  1  0.4468985    
## age:co2:salinity                 0.0715  1  0.7891903    
## age:co2:elevation                5.4149  1  0.0199660 *  
## age:salinity:elevation           1.2144  1  0.2704696    
## location:co2:salinity            0.0891  1  0.7652656    
## location:co2:elevation           0.0846  1  0.7711070    
## location:salinity:elevation      2.1855  1  0.1393118    
## co2:salinity:elevation           0.1411  1  0.7071706    
## age:location:co2:salinity        0.1081  1  0.7423733    
## age:location:co2:elevation       0.5091  1  0.4755231    
## age:location:salinity:elevation  0.0940  1  0.7591596    
## age:co2:salinity:elevation       0.0103  1  0.9190293    
## location:co2:salinity:elevation  0.0878  1  0.7669450    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

rs_mod3 <- lmer(log(rs) ~ weight_init + date_planted_grp + origin_lab +
                  (age + location + co2 + salinity + elevation)^3 + I(elevation^2) +
                  + (1|site_frame) + (1|genotype), data = bg_nocomp_rs)

step(rs_mod3, keep = c("co2", "salinity", "elevation", "age", "location"))

## Backward reduced random-effect table:
## 
##                  Eliminated npar  logLik    AIC     LRT Df Pr(>Chisq)    
## <none>                        33 -23.116 112.23                          
## (1 | site_frame)          1   32 -23.367 110.73  0.5025  1     0.4784    
## (1 | genotype)            0   31 -31.591 125.18 16.4477  1  5.001e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                             Eliminated  Sum Sq Mean Sq NumDF   DenDF F value
## location:co2:elevation               1 0.00470 0.00470     1 182.373  0.0883
## location:co2:salinity                2 0.00699 0.00699     1 176.065  0.1322
## co2:salinity:elevation               3 0.00646 0.00646     1 190.483  0.1227
## age:co2:salinity                     4 0.00769 0.00769     1 177.362  0.1470
## origin_lab                           5 0.01622 0.01622     1  34.870  0.3115
## weight_init                          6 0.02904 0.02904     1 194.177  0.5569
## age:location:co2                     7 0.03373 0.03373     1 186.721  0.6473
## location:co2                         8 0.00016 0.00016     1 187.020  0.0030
## co2:salinity                         9 0.03929 0.03929     1 181.555  0.7608
## age:location:elevation              10 0.04495 0.04495     1 186.711  0.8727
## age:salinity:elevation              11 0.05752 0.05752     1 187.427  1.1129
## date_planted_grp                    12 0.05223 0.05223     1 193.891  1.0090
## age:location:salinity               13 0.08469 0.08469     1 191.077  1.6389
## age:location                        14 0.00834 0.00834     1  25.358  0.1610
## age:salinity                        15 0.03898 0.03898     1 193.725  0.7523
## location:salinity:elevation         16 0.09355 0.09355     1 193.763  1.8097
## location:salinity                   17 0.11842 0.11842     1 193.828  2.2880
## I(elevation^2)                      18 0.17629 0.17629     1 196.678  3.3896
## location:elevation                   0 0.26582 0.26582     1 194.340  5.0620
## salinity:elevation                   0 0.71603 0.71603     1 194.686 13.6353
## age:co2:elevation                    0 0.28637 0.28637     1 199.725  5.4533
##                                Pr(>F)    
## location:co2:elevation      0.7666359    
## location:co2:salinity       0.7165510    
## co2:salinity:elevation      0.7265044    
## age:co2:salinity            0.7018869    
## origin_lab                  0.5803156    
## weight_init                 0.4564266    
## age:location:co2            0.4221042    
## location:co2                0.9560358    
## co2:salinity                0.3842361    
## age:location:elevation      0.3514056    
## age:salinity:elevation      0.2928045    
## date_planted_grp            0.3164039    
## age:location:salinity       0.2020318    
## age:location                0.6915659    
## age:salinity                0.3868097    
## location:salinity:elevation 0.1801213    
## location:salinity           0.1320070    
## I(elevation^2)              0.0671155 .  
## location:elevation          0.0255772 *  
## salinity:elevation          0.0002882 ***
## age:co2:elevation           0.0205239 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## log(rs) ~ age + location + co2 + salinity + elevation + (1 | 
##     genotype) + age:co2 + age:elevation + location:elevation + 
##     co2:elevation + salinity:elevation + age:co2:elevation

# Create the stepwise model
rs_mod4 <- lmer(log(rs) ~ location*elevation + salinity*elevation + age*co2*elevation+
                 (1|genotype), data = bg_nocomp_rs)

# Diagnostic plots
plot_model(rs_mod4, type = "diag")

## [[1]]

## `geom_smooth()` using formula 'y ~ x'

## 
## [[2]]
## [[2]]$genotype

## `geom_smooth()` using formula 'y ~ x'

## 
## 
## [[3]]

## 
## [[4]]

## `geom_smooth()` using formula 'y ~ x'

# See if there are random effects models that are an improvement of this model
rs_mod5 <- update(rs_mod4, .~.-(1|genotype) + (co2*elevation|genotype))

## boundary (singular) fit: see ?isSingular

rs_mod6 <- update(rs_mod4, .~.-(1|genotype) + (co2*salinity|genotype))

## boundary (singular) fit: see ?isSingular

## Warning: Model failed to converge with 1 negative eigenvalue: -8.2e-03

rs_mod7 <- update(rs_mod4, .~.-(1|genotype) + (salinity*elevation|genotype)) # fits
rs_mod8 <- update(rs_mod4, .~.-(1|genotype) + (co2+elevation|genotype))

## boundary (singular) fit: see ?isSingular

## Warning: Model failed to converge with 1 negative eigenvalue: -9.4e-05

rs_mod9 <- update(rs_mod4, .~.-(1|genotype) + (co2+salinity|genotype))
rs_mod10 <- update(rs_mod4, .~.-(1|genotype) + (salinity+elevation|genotype)) # fits
rs_mod11 <- update(rs_mod4, .~.-(1|genotype) + (co2|genotype)) 

## boundary (singular) fit: see ?isSingular

rs_mod12 <- update(rs_mod4, .~.-(1|genotype) + (salinity|genotype)) # fits 
rs_mod13 <- update(rs_mod4, .~.-(1|genotype) + (elevation|genotype)) # fits

anova(rs_mod4, rs_mod7, refit = F) # ns

## Data: bg_nocomp_rs
## Models:
## rs_mod4: log(rs) ~ location * elevation + salinity * elevation + age * 
## rs_mod4:     co2 * elevation + (1 | genotype)
## rs_mod7: log(rs) ~ location + elevation + salinity + age + co2 + (salinity * 
## rs_mod7:     elevation | genotype) + location:elevation + elevation:salinity + 
## rs_mod7:     age:co2 + elevation:age + elevation:co2 + elevation:age:co2
##         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## rs_mod4   14 57.793 105.68 -14.897   29.793                     
## rs_mod7   23 67.283 145.96 -10.642   21.283 8.5103  9     0.4836

anova(rs_mod4, rs_mod10, refit = F) # ns

## Data: bg_nocomp_rs
## Models:
## rs_mod4: log(rs) ~ location * elevation + salinity * elevation + age * 
## rs_mod4:     co2 * elevation + (1 | genotype)
## rs_mod10: log(rs) ~ location + elevation + salinity + age + co2 + (salinity + 
## rs_mod10:     elevation | genotype) + location:elevation + elevation:salinity + 
## rs_mod10:     age:co2 + elevation:age + elevation:co2 + elevation:age:co2
##          npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## rs_mod4    14 57.793 105.68 -14.897   29.793                    
## rs_mod10   19 62.433 127.42 -12.217   24.433  5.36  5     0.3735

anova(rs_mod4, rs_mod12, refit = F) # ns

## Data: bg_nocomp_rs
## Models:
## rs_mod4: log(rs) ~ location * elevation + salinity * elevation + age * 
## rs_mod4:     co2 * elevation + (1 | genotype)
## rs_mod12: log(rs) ~ location + elevation + salinity + age + co2 + (salinity | 
## rs_mod12:     genotype) + location:elevation + elevation:salinity + age:co2 + 
## rs_mod12:     elevation:age + elevation:co2 + elevation:age:co2
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## rs_mod4    14 57.793 105.68 -14.897   29.793                     
## rs_mod12   16 58.109 112.84 -13.055   26.109 3.6842  2     0.1585

anova(rs_mod4, rs_mod13, refit = F) # ns

## Data: bg_nocomp_rs
## Models:
## rs_mod4: log(rs) ~ location * elevation + salinity * elevation + age * 
## rs_mod4:     co2 * elevation + (1 | genotype)
## rs_mod13: log(rs) ~ location + elevation + salinity + age + co2 + (elevation | 
## rs_mod13:     genotype) + location:elevation + elevation:salinity + age:co2 + 
## rs_mod13:     elevation:age + elevation:co2 + elevation:age:co2
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## rs_mod4    14 57.793 105.68 -14.897   29.793                     
## rs_mod13   16 60.732 115.46 -14.366   28.732 1.0617  2     0.5881

# None of the random effects add significantly

# Create caterpillar graph
REsim(rs_mod4) %>%
  mutate(groupID = fct_reorder(groupID, mean)) %>% 
  ggplot(aes(x = groupID)) +
  geom_pointrange(aes(ymin = mean - 2*sd, ymax = mean + 2*sd, xmin = groupID, xmax = groupID, y = mean, color = groupID)) +
  coord_flip() +
  ylab("Random effect") +
  xlab("Genotype") +
  geom_hline(aes(yintercept = 0), linetype = "dashed") +
  theme_bw() + scale_x_discrete(position = "top") +
  theme(legend.position = "none") +
  scale_color_manual(values = col_vector)-> d

# Also fit a null model
rs_mod_null <- lm(log(rs) ~ elevation*salinity + co2*elevation, data = bg_nocomp_rs)

# Compare R2 values
rs_null_r2 <- r.squaredGLMM(rs_mod_null)[2]
rs_r2 <- r.squaredGLMM(rs_mod4)[2]

# Make plots of fixed effects
Anova(rs_mod4, type = 3)

## Analysis of Deviance Table (Type III Wald chisquare tests)
## 
## Response: log(rs)
##                      Chisq Df Pr(>Chisq)    
## (Intercept)         7.4957  1   0.006185 ** 
## location            2.1537  1   0.142231    
## elevation          56.3481  1  6.071e-14 ***
## salinity           18.6729  1  1.552e-05 ***
## age                 0.4679  1   0.493940    
## co2                 0.4316  1   0.511207    
## location:elevation  5.0620  1   0.024456 *  
## elevation:salinity 13.6353  1   0.000222 ***
## age:co2             4.9595  1   0.025948 *  
## elevation:age       0.3915  1   0.531489    
## elevation:co2       2.0213  1   0.155106    
## elevation:age:co2   5.4533  1   0.019531 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Need elevation:age:co2, elevation:salinity, and location:elevation

# Extract predicted means for model
rs_means <- emmeans::emmeans(rs_mod4, ~co2:age:elevation,
                             at= list(elevation = seq(0.15, 0.55, 0.01)), type = "response")

rs_means_out <- summary(rs_means)

rs_means_out %>% 
  mutate(age = case_when(age == "modern" ~ "Descendant cohort (2000-2020)",
                         T ~ "Ancestral cohort (1900-1960)")) -> rs_means_out

bg_nocomp_rs_plot <- bg_nocomp_plot %>% 
  mutate(rs = total_bg / agb_scam) %>% 
  filter(rs < 5)

rs_means_out %>% 
  ggplot(aes(x = elevation, y = response, color = co2)) + 
  geom_line(size = 1.5) +
  facet_wrap(~age) +
  geom_point(data = bg_nocomp_rs_plot, aes(x = elevation, y = rs), alpha = 0.5) +
  theme_bw() +
  scale_color_manual(values = c("gray47", "#66a61e")) +
  ylab("Root-to-shoot ratio") +
  xlab("Elevation (m NAVD88)") +
  guides(color=guide_legend(title=expression(CO[2]))) +
  theme(legend.position = "top") -> e

rs_means2 <- emmeans(rs_mod4, ~salinity:elevation,
                    at= list(elevation = seq(0.15, 0.55, 0.01)), type = "response")

rs_means_out2 <- summary(rs_means2)

rs_means_out2 %>% 
  ggplot(aes(x = elevation, y = response, color = salinity)) +
  geom_line(size = 1.5) +
  geom_point(data = bg_nocomp_rs_plot, aes(x = elevation, y = rs), alpha = 0.5) +
  theme_bw() + scale_color_manual(values = c("orange", "gray47")) +
  ylab("Root-to-shoot ratio") +
  xlab("Elevation (m NAVD88)") +
  guides(color=guide_legend("")) +
  theme(legend.position = "top") -> f

rs_means3 <- emmeans(rs_mod4, ~location:elevation,
                    at= list(elevation = seq(0.15, 0.55, 0.01)), type = "response")

rs_means_out3 <- summary(rs_means3)

rs_means_out3 %>% 
  mutate(location = case_when(location == "corn" ~ "Corn Island",
                              T ~ "Kirkpatrick Marsh")) -> rs_means_out3

bg_nocomp_rs_plot %>% 
  mutate(location = case_when(location == "corn" ~ "Corn Island",
                              T ~ "Kirkpatrick Marsh")) -> bg_nocomp_rs_plot

rs_means_out3 %>% 
  ggplot(aes(x = elevation, y = response, color = location)) +
  geom_line(size = 1.5) +
  geom_point(data = bg_nocomp_rs_plot, aes(x = elevation, y = rs), alpha = 0.5) +
  theme_bw() + scale_color_manual(values = c("#1b9e77", "#7570b3")) +
  ylab("Root-to-shoot ratio") +
  xlab("Elevation (m NAVD88)") +
  guides(color=guide_legend("")) +
  theme(legend.position = "top") -> g

cowplot::plot_grid(f, g, labels = c("b", "c")) -> fg
cowplot::plot_grid(e, fg, ncol = 1, labels = c("a", "", "")) -> fge
cowplot::plot_grid(d, labels = "d") -> d
fge + d + plot_layout(widths = c(4,1), guides = "collect") -> rs_plot

Root distribution parameter

# Most complex model
beta_mod <- lmer(beta ~ weight_init + date_planted_grp + origin_lab +
                 (age + location + co2 + salinity + elevation_sc)^5 + I(elevation_sc^2) +
                 (1|site_frame) + (1|genotype), data = bg_nocomp_rs)

# Do stepwise model selection
# step(beta_mod) # too complex

VarCorr(beta_mod)

##  Groups     Name        Std.Dev. 
##  genotype   (Intercept) 0.0119866
##  site_frame (Intercept) 0.0032349
##  Residual               0.0266574

Anova(beta_mod)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: beta
##                                           Chisq Df Pr(>Chisq)    
## weight_init                              1.8577  1    0.17289    
## date_planted_grp                         0.2637  1    0.60761    
## origin_lab                               0.8143  1    0.36686    
## age                                      3.6133  1    0.05732 .  
## location                                 0.3393  1    0.56025    
## co2                                      1.1127  1    0.29150    
## salinity                                 0.5640  1    0.45264    
## elevation_sc                           204.6761  1    < 2e-16 ***
## I(elevation_sc^2)                        2.4456  1    0.11785    
## age:location                             0.0324  1    0.85722    
## age:co2                                  0.2908  1    0.58973    
## age:salinity                             1.0714  1    0.30062    
## age:elevation_sc                         1.5983  1    0.20614    
## location:co2                             0.0435  1    0.83487    
## location:salinity                        4.0349  1    0.04457 *  
## location:elevation_sc                    4.2626  1    0.03896 *  
## co2:salinity                             0.5169  1    0.47215    
## co2:elevation_sc                         4.2644  1    0.03892 *  
## salinity:elevation_sc                    5.9665  1    0.01458 *  
## age:location:co2                         1.1063  1    0.29288    
## age:location:salinity                    0.6866  1    0.40731    
## age:location:elevation_sc                0.3536  1    0.55207    
## age:co2:salinity                         0.0112  1    0.91580    
## age:co2:elevation_sc                     2.8633  1    0.09062 .  
## age:salinity:elevation_sc                1.2128  1    0.27079    
## location:co2:salinity                    0.1042  1    0.74688    
## location:co2:elevation_sc                0.4084  1    0.52277    
## location:salinity:elevation_sc           1.0313  1    0.30984    
## co2:salinity:elevation_sc                2.7875  1    0.09500 .  
## age:location:co2:salinity                0.0041  1    0.94911    
## age:location:co2:elevation_sc            0.6246  1    0.42933    
## age:location:salinity:elevation_sc       1.4932  1    0.22172    
## age:co2:salinity:elevation_sc            0.1184  1    0.73074    
## location:co2:salinity:elevation_sc       0.1714  1    0.67886    
## age:location:co2:salinity:elevation_sc   0.1397  1    0.70859    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# site_frame variance is low and no significant 5-way and 4-way interactions
beta_mod2 <- update(beta_mod, .~.-(1|site_frame))
# step(beta_mod2)

beta_mod3 <- lmer(beta ~ weight_init + date_planted_grp + origin_lab +
                 (age + location + co2 + salinity + elevation_sc)^4 + I(elevation_sc^2) +
                 (1|genotype), data = bg_nocomp_rs)

step(beta_mod3, keep = c("co2", "salinity", "elevation_sc", "age", "location"))

## Backward reduced random-effect table:
## 
##                Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)   
## <none>                      37 355.15 -636.30                        
## (1 | genotype)          0   36 350.49 -628.98 9.3233  1   0.002263 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                                    Eliminated    Sum Sq   Mean Sq NumDF   DenDF
## age:location:co2:salinity                   1 0.0000035 0.0000035     1 169.793
## age:co2:salinity:elevation_sc               2 0.0000897 0.0000897     1 182.482
## age:co2:salinity                            3 0.0000147 0.0000147     1 169.454
## location:co2:salinity:elevation_sc          4 0.0001251 0.0001251     1 190.543
## location:co2:salinity                       5 0.0000678 0.0000678     1 171.902
## date_planted_grp                            6 0.0001786 0.0001786     1 183.350
## origin_lab                                  7 0.0005246 0.0005246     1  33.553
## age:location:co2:elevation_sc               8 0.0004973 0.0004973     1 182.509
## location:co2:elevation_sc                   9 0.0002324 0.0002324     1 182.652
## age:location:co2                           10 0.0008218 0.0008218     1 182.753
## location:co2                               11 0.0000119 0.0000119     1 183.807
## age:location:salinity:elevation_sc         12 0.0008337 0.0008337     1 183.998
## age:location:elevation_sc                  13 0.0001855 0.0001855     1 183.555
## age:location:salinity                      14 0.0004579 0.0004579     1 184.733
## age:location                               15 0.0000516 0.0000516     1  26.153
## location:salinity:elevation_sc             16 0.0007878 0.0007878     1 188.462
## age:salinity:elevation_sc                  17 0.0008435 0.0008435     1 186.440
## age:salinity                               18 0.0005978 0.0005978     1 189.759
## weight_init                                19 0.0007834 0.0007834     1 201.412
## co2:salinity:elevation_sc                  20 0.0013889 0.0013889     1 200.890
## co2:salinity                               21 0.0006082 0.0006082     1 186.612
## I(elevation_sc^2)                          22 0.0015584 0.0015584     1 194.340
## age:co2:elevation_sc                       23 0.0025997 0.0025997     1 197.628
## age:co2                                    24 0.0002661 0.0002661     1 197.664
## age:elevation_sc                           25 0.0013573 0.0013573     1 192.175
## co2:elevation_sc                           26 0.0020924 0.0020924     1 199.335
## location:salinity                          27 0.0025404 0.0025404     1 196.318
## location:elevation_sc                      28 0.0026331 0.0026331     1 195.520
## age                                         0 0.0026093 0.0026093     1  26.707
## location                                    0 0.0005226 0.0005226     1  27.041
## co2                                         0 0.0008293 0.0008293     1 202.669
## salinity:elevation_sc                       0 0.0039719 0.0039719     1 198.555
##                                    F value  Pr(>F)  
## age:location:co2:salinity           0.0049 0.94400  
## age:co2:salinity:elevation_sc       0.1263 0.72267  
## age:co2:salinity                    0.0208 0.88561  
## location:co2:salinity:elevation_sc  0.1783 0.67336  
## location:co2:salinity               0.0972 0.75558  
## date_planted_grp                    0.2575 0.61242  
## origin_lab                          0.7591 0.38981  
## age:location:co2:elevation_sc       0.7203 0.39716  
## location:co2:elevation_sc           0.3374 0.56204  
## age:location:co2                    1.1999 0.27478  
## location:co2                        0.0172 0.89568  
## age:location:salinity:elevation_sc  1.2178 0.27123  
## age:location:elevation_sc           0.2713 0.60309  
## age:location:salinity               0.6732 0.41298  
## age:location                        0.0761 0.78482  
## location:salinity:elevation_sc      1.1625 0.28232  
## age:salinity:elevation_sc           1.2465 0.26566  
## age:salinity                        0.8816 0.34895  
## weight_init                         1.1559 0.28360  
## co2:salinity:elevation_sc           2.0460 0.15416  
## co2:salinity                        0.8908 0.34648  
## I(elevation_sc^2)                   2.2878 0.13202  
## age:co2:elevation_sc                3.7853 0.05312 .
## age:co2                             0.3853 0.53548  
## age:elevation_sc                    1.9713 0.16192  
## co2:elevation_sc                    3.0235 0.08361 .
## location:salinity                   3.6577 0.05727 .
## location:elevation_sc               3.7608 0.05391 .
## age                                 3.6662 0.06629 .
## location                            0.7343 0.39902  
## co2                                 1.1652 0.28167  
## salinity:elevation_sc               5.5807 0.01913 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## beta ~ age + location + co2 + salinity + elevation_sc + (1 | 
##     genotype) + salinity:elevation_sc

# Create stepwise model
beta_mod4 <- lmer(beta ~ age + location + co2 + salinity*elevation_sc +
                    (1 | genotype), data = bg_nocomp_rs)

# Check diagnostics
plot_model(beta_mod4, type = "diag") 

## [[1]]

## `geom_smooth()` using formula 'y ~ x'

## 
## [[2]]
## [[2]]$genotype

## `geom_smooth()` using formula 'y ~ x'

## 
## 
## [[3]]

## 
## [[4]]

## `geom_smooth()` using formula 'y ~ x'

# Some evidence of non-normality, but not horrible

# See if there are significant random intercepts
beta_mod5 <- lmer(beta ~ age + location + co2 + salinity*elevation_sc +
                    (1 + salinity*elevation_sc| genotype), data = bg_nocomp_rs) # fits
beta_mod6 <- lmer(beta ~ age + location + co2 + salinity*elevation_sc +
                    (1 + salinity+elevation_sc| genotype), data = bg_nocomp_rs) # fits
beta_mod7 <- lmer(beta ~ age + location + co2 + salinity*elevation_sc +
                    (1 + salinity+co2| genotype), data = bg_nocomp_rs)

## boundary (singular) fit: see ?isSingular

beta_mod8 <- lmer(beta ~ age + location + co2 + salinity*elevation_sc +
                    (1 + elevation_sc+co2| genotype), data = bg_nocomp_rs)

## boundary (singular) fit: see ?isSingular

beta_mod9 <- lmer(beta ~ age + location + co2 + salinity*elevation_sc +
                    (1 + co2| genotype), data = bg_nocomp_rs)

## boundary (singular) fit: see ?isSingular

beta_mod10 <- lmer(beta ~ age + location + co2 + salinity*elevation_sc +
                    (1 + salinity| genotype), data = bg_nocomp_rs) # fits
beta_mod11 <- lmer(beta ~ age + location + co2 + salinity*elevation_sc +
                    (1 + elevation_sc| genotype), data = bg_nocomp_rs) # fits

anova(beta_mod5, beta_mod4) # *** best model

## refitting model(s) with ML (instead of REML)

## Data: bg_nocomp_rs
## Models:
## beta_mod4: beta ~ age + location + co2 + salinity * elevation_sc + (1 | 
## beta_mod4:     genotype)
## beta_mod5: beta ~ age + location + co2 + salinity * elevation_sc + (1 + 
## beta_mod5:     salinity * elevation_sc | genotype)
##           npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)    
## beta_mod4    9 -940.68 -910.02 479.34  -958.68                         
## beta_mod5   18 -958.06 -896.73 497.03  -994.06 35.371  9  5.124e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(beta_mod6, beta_mod5) # ***

## refitting model(s) with ML (instead of REML)

## Data: bg_nocomp_rs
## Models:
## beta_mod6: beta ~ age + location + co2 + salinity * elevation_sc + (1 + 
## beta_mod6:     salinity + elevation_sc | genotype)
## beta_mod5: beta ~ age + location + co2 + salinity * elevation_sc + (1 + 
## beta_mod5:     salinity * elevation_sc | genotype)
##           npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)  
## beta_mod6   14 -955.93 -908.23 491.97  -983.93                       
## beta_mod5   18 -958.06 -896.73 497.03  -994.06 10.122  4    0.03843 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(beta_mod10, beta_mod4) # *

## refitting model(s) with ML (instead of REML)

## Data: bg_nocomp_rs
## Models:
## beta_mod4: beta ~ age + location + co2 + salinity * elevation_sc + (1 | 
## beta_mod4:     genotype)
## beta_mod10: beta ~ age + location + co2 + salinity * elevation_sc + (1 + 
## beta_mod10:     salinity | genotype)
##            npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)  
## beta_mod4     9 -940.68 -910.02 479.34  -958.68                       
## beta_mod10   11 -943.92 -906.44 482.96  -965.92 7.2339  2    0.02686 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(beta_mod11, beta_mod4) # **

## refitting model(s) with ML (instead of REML)

## Data: bg_nocomp_rs
## Models:
## beta_mod4: beta ~ age + location + co2 + salinity * elevation_sc + (1 | 
## beta_mod4:     genotype)
## beta_mod11: beta ~ age + location + co2 + salinity * elevation_sc + (1 + 
## beta_mod11:     elevation_sc | genotype)
##            npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)   
## beta_mod4     9 -940.68 -910.02 479.34  -958.68                        
## beta_mod11   11 -946.63 -909.15 484.31  -968.63 9.9407  2   0.006941 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Best model has salinity*elevation random slope
beta_finalmod <-lmer(beta ~ age + location + co2 + salinity*elevation_sc +
                    (1 + salinity*elevation_sc| genotype), data = bg_nocomp_rs)

# Look at random effect
out <- plot_model(beta_finalmod, type = "pred", pred.type = "re", terms = c("elevation_sc[all]", "genotype","salinity")) +
  scale_color_manual(values = rainbow(31)) +
  ggtitle("") +
  theme_bw() +
  theme(legend.position = "none") +
  ylab("Root distribution parameter") +
  xlab("Elevation (m NAVD88)") 

## Warning in RColorBrewer::brewer.pal(n, pal): n too large, allowed maximum for palette Set1 is 9
## Returning the palette you asked for with that many colors

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

# Extract values from plot_model so we can customize
tibble(elevation_sc = out$data$x,
       root_param = out$data$predicted,
       genotype = out$data$group,
       salinity = out$data$facet) %>% 
  mutate(elevation = elevation_sc*scaling + centering) -> beta_GxE_data

beta_GxE_data %>% 
  filter(elevation %in% min(elevation)) -> mins

beta_GxE_data %>% 
  filter(elevation %in% max(elevation)) -> maxs

# Make random effects plot
beta_GxE_data %>% 
  ggplot(aes(x = elevation, y = root_param, color = genotype)) +
  geom_line(size = 1, alpha = 1) +
  facet_wrap(~salinity) +
  scale_color_manual(values = col_vector) +
  theme_bw() +
  theme(legend.position = "none") +
  ylab("Root distribution parameter") +
  xlab("Elevation (m NAVD88)") +
  geom_point(data = mins, aes(x = elevation, y = root_param, color = genotype), size = 2.7) +
  geom_point(data = maxs, aes(x = elevation, y = root_param, color = genotype), size = 2.7) -> h
 

# Make fixed effect plot
beta_means <- emmeans(beta_finalmod, ~salinity:elevation_sc,
                    at= list(elevation_sc = seq(-1.72, 1.54, 0.1)), type = "response")

beta_means_out <- summary(beta_means)

beta_means_out %>% 
  mutate(elevation = elevation_sc*scaling + centering) %>% 
  ggplot(aes(x = elevation, y = emmean, color = salinity)) +
  geom_line(size = 1.5) +
  geom_point(data = bg_nocomp_rs_plot, aes(x = elevation, y = beta), alpha = 0.5) +
  theme_bw() + scale_color_manual(values = c("orange", "gray47")) +
  ylab("Root distribution parameter") +
  xlab("Elevation (m NAVD88)") +
  guides(color=guide_legend("")) +
  theme(legend.position = "top") +
  annotate("text", x = c(0.22, 0.22), y = c(0.95, 0.77), label = c("deeper rooting", "shallower rooting"),color = "gray47", fontface = 3) -> i


beta_finalmod <-lmer(beta ~ age + location + co2 + salinity*elevation_sc +
                    (1 + salinity*elevation_sc| genotype), data = bg_nocomp_rs_plot)

my_labeller <- as_labeller(function(x){
  return(paste0(x))
})

# Make fixed effect plot for age overall
emmip(beta_finalmod, location ~ age, CI = T, weights = "proportional") +
  geom_jitter(data = bg_nocomp_rs_plot, aes(x = age, y = beta, color = location),
              height = 0, width = 0.1, alpha = 0.3, shape = 1) +
  ylab("Root distribution parameter") +
  scale_x_discrete("Age cohort",
                   labels = c("Descendant cohort (2000-2020)" = "Descendant (2000-2020)",
                              "Ancestral cohort (1900-1960)" = "Ancestral (1900-1960)")) +
  scale_color_manual(values = c("#1b9e77", "#7570b3")) + theme_bw() +
  theme(legend.position = "top") + guides(color=guide_legend("")) +
  scale_size_manual(values=c(2,2)) -> j

(i + j) / h + plot_annotation(tag_levels = 'a') -> beta_plot

# Compare to a null model
beta_null_mod <-lm(beta ~ co2 + salinity*elevation_sc, data = bg_nocomp_rs_plot)
beta_null_r2 <- r.squaredGLMM(beta_null_mod)[2]
beta_r2 <- r.squaredGLMM(beta_finalmod)[2]

Average stem height

height_mod <- lmer(height_scam_tot ~ weight_init_sc + date_planted_grp + origin_lab +
                  (age + location + co2 + salinity + elevation_sc)^5 + I(elevation_sc^2) +
                  (1|site_frame) + (1|genotype), data = bg_nocomp)

plot(height_mod)# Major outlier

bg_nocomp %>% filter(pot_no == 1830)

##   pot_no    beta bg_scam10 bg_sppa10 total_bg site_frame level position
## 1   1830 0.93414     0.637        NA    1.239    gcrew_5     1        2
##        co2 salinity comp elevation genotype location depth_seed         cohort
## 1 elevated     salt    0 0.5263333      ca7     corn      20.75 corn_ancestral
##   origin_lab date_cloned_grp date_planted_grp weight_init width_init nodes_init
## 1        utk               1                2       0.948       1.69          1
##   agb_scam agb_sppa dens_scam_live dens_scam_dead height_sppa height_scam_tot
## 1    0.721        0              7              1          NA            11.5
##   height_scam_grn width_scam_mid width_scam_bas prop_twist soil_to_pot
## 1            11.5            1.7            1.7          0           6
##            id extinct       age elevation_sc weight_init_sc       rs
## 1 gcrew_5_1_2       0 ancestral     1.395559     -0.3406356 1.718447

# This is the outlier and note on data sheet says that this pot had significant
# herbivory which would affect the average height

bg_nocomp_height <- bg_nocomp %>% 
  filter(pot_no != 1830)

height_mod2 <- lmer(height_scam_tot ~ weight_init_sc + date_planted_grp + origin_lab +
                  (age + location + co2 + salinity + elevation_sc)^5 + I(elevation_sc^2) +
                  (1|site_frame) + (1|genotype), data = bg_nocomp_height)

# Stepwise selection
#step(height_mod2) # Too complex

VarCorr(height_mod2)

##  Groups     Name        Std.Dev.
##  genotype   (Intercept) 3.57274 
##  site_frame (Intercept) 0.34105 
##  Residual               4.62926

Anova(height_mod2)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: height_scam_tot
##                                           Chisq Df Pr(>Chisq)    
## weight_init_sc                           2.5785  1    0.10832    
## date_planted_grp                         2.5936  1    0.10730    
## origin_lab                               0.1670  1    0.68278    
## age                                      0.1072  1    0.74338    
## location                                 0.0539  1    0.81645    
## co2                                      0.0164  1    0.89818    
## salinity                                 0.1062  1    0.74455    
## elevation_sc                           147.9034  1    < 2e-16 ***
## I(elevation_sc^2)                        0.2943  1    0.58747    
## age:location                             0.5482  1    0.45904    
## age:co2                                  0.5388  1    0.46291    
## age:salinity                             0.5105  1    0.47494    
## age:elevation_sc                         0.9837  1    0.32128    
## location:co2                             0.7792  1    0.37739    
## location:salinity                        3.8988  1    0.04832 *  
## location:elevation_sc                    2.4535  1    0.11726    
## co2:salinity                             0.3058  1    0.58028    
## co2:elevation_sc                         0.8110  1    0.36781    
## salinity:elevation_sc                    0.5857  1    0.44408    
## age:location:co2                         0.8210  1    0.36488    
## age:location:salinity                    0.1151  1    0.73437    
## age:location:elevation_sc                1.4778  1    0.22412    
## age:co2:salinity                         0.1942  1    0.65940    
## age:co2:elevation_sc                     1.0207  1    0.31235    
## age:salinity:elevation_sc                1.1686  1    0.27968    
## location:co2:salinity                    0.3695  1    0.54330    
## location:co2:elevation_sc                0.2627  1    0.60830    
## location:salinity:elevation_sc           0.1723  1    0.67808    
## co2:salinity:elevation_sc                0.7665  1    0.38131    
## age:location:co2:salinity                1.3312  1    0.24859    
## age:location:co2:elevation_sc            1.7941  1    0.18043    
## age:location:salinity:elevation_sc       0.4387  1    0.50776    
## age:co2:salinity:elevation_sc            1.1499  1    0.28358    
## location:co2:salinity:elevation_sc       0.3212  1    0.57086    
## age:location:co2:salinity:elevation_sc   0.5296  1    0.46677    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Drop 5-way interaction
height_mod3 <- lmer(height_scam_tot ~ weight_init_sc + date_planted_grp + origin_lab +
                  (age + location + co2 + salinity + elevation_sc)^4 + I(elevation_sc^2) +
                  (1|site_frame) + (1|genotype), data = bg_nocomp_height)

#step(height_mod3)
VarCorr(height_mod3)

##  Groups     Name        Std.Dev.
##  genotype   (Intercept) 3.60979 
##  site_frame (Intercept) 0.35074 
##  Residual               4.61663

height_mod4 <- lmer(height_scam_tot ~ weight_init_sc + date_planted_grp + origin_lab +
                  (age + location + co2 + salinity + elevation_sc)^4 +
                    +I(elevation_sc^2) + (1|genotype), data = bg_nocomp_height)

step(height_mod4, keep = c("co2", "salinity", "age", "location", "elevation_sc"))

## Backward reduced random-effect table:
## 
##                Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                      37 -638.59 1351.2                         
## (1 | genotype)          0   36 -662.91 1397.8 48.642  1  3.072e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                                    Eliminated Sum Sq Mean Sq NumDF   DenDF
## origin_lab                                  1    3.4     3.4     1  42.738
## I(elevation_sc^2)                           2    5.6     5.6     1 177.115
## location:co2:salinity:elevation_sc          3    6.9     6.9     1 189.353
## age:location:salinity:elevation_sc          4    7.2     7.2     1 176.515
## location:salinity:elevation_sc              5    3.4     3.4     1 177.694
## age:co2:salinity:elevation_sc               6   24.0    24.0     1 186.491
## co2:salinity:elevation_sc                   7   14.2    14.2     1 186.332
## age:salinity:elevation_sc                   8   22.0    22.0     1 179.905
## salinity:elevation_sc                       9   11.2    11.2     1 181.159
## age:location:co2:salinity                  10   25.6    25.6     1 180.483
## age:location:salinity                      11    2.2     2.2     1 183.906
## age:co2:salinity                           12    6.1     6.1     1 182.057
## location:co2:salinity                      13    9.9     9.9     1 183.148
## co2:salinity                               14    5.8     5.8     1 183.800
## age:salinity                               15   12.1    12.1     1 188.737
## age:location:co2:elevation_sc              16   26.5    26.5     1 190.536
## location:co2:elevation_sc                  17    4.1     4.1     1 192.002
## age:location:co2                           18   17.1    17.1     1 191.656
## location:co2                               19   19.5    19.5     1 193.314
## age:co2:elevation_sc                       20   25.1    25.1     1 195.232
## age:co2                                    21    8.4     8.4     1 196.154
## co2:elevation_sc                           22   17.6    17.6     1 197.555
## age:location:elevation_sc                  23   38.3    38.3     1 194.616
## age:location                               24   10.7    10.7     1  28.090
## age:elevation_sc                           25   19.2    19.2     1 195.608
## location:elevation_sc                      26   31.7    31.7     1 196.096
## weight_init_sc                             27   40.1    40.1     1 206.394
## location:salinity                          28   70.2    70.2     1 200.231
## date_planted_grp                            0  145.9   145.9     1 202.310
## age                                         0    1.6     1.6     1  29.057
## location                                    0    2.1     2.1     1  29.310
## co2                                         0    2.1     2.1     1 203.494
## salinity                                    0   18.6    18.6     1 202.038
## elevation_sc                                0 3429.1  3429.1     1 198.957
##                                     F value    Pr(>F)    
## origin_lab                           0.1588  0.692275    
## I(elevation_sc^2)                    0.2616  0.609671    
## location:co2:salinity:elevation_sc   0.3230  0.570484    
## age:location:salinity:elevation_sc   0.3405  0.560301    
## location:salinity:elevation_sc       0.1619  0.687876    
## age:co2:salinity:elevation_sc        1.1398  0.287079    
## co2:salinity:elevation_sc            0.6754  0.412225    
## age:salinity:elevation_sc            1.0475  0.307451    
## salinity:elevation_sc                0.5337  0.466014    
## age:location:co2:salinity            1.2224  0.270353    
## age:location:salinity                0.1063  0.744802    
## age:co2:salinity                     0.2907  0.590431    
## location:co2:salinity                0.4762  0.491030    
## co2:salinity                         0.2781  0.598563    
## age:salinity                         0.5844  0.445557    
## age:location:co2:elevation_sc        1.2846  0.258471    
## location:co2:elevation_sc            0.1965  0.658048    
## age:location:co2                     0.8318  0.362889    
## location:co2                         0.9496  0.331046    
## age:co2:elevation_sc                 1.2178  0.271159    
## age:co2                              0.4051  0.525202    
## co2:elevation_sc                     0.8533  0.356736    
## age:location:elevation_sc            1.8618  0.173992    
## age:location                         0.5202  0.476699    
## age:elevation_sc                     0.9287  0.336387    
## location:elevation_sc                1.5365  0.216616    
## weight_init_sc                       1.9409  0.165069    
## location:salinity                    3.3881  0.067148 .  
## date_planted_grp                     6.9633  0.008967 ** 
## age                                  0.0758  0.784991    
## location                             0.0986  0.755695    
## co2                                  0.1012  0.750707    
## salinity                             0.8867  0.347500    
## elevation_sc                       163.6906 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + 
##     elevation_sc + (1 | genotype)

# Fit proposed model from step procedure
height_mod5 <- lmer(height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + elevation_sc +
                      (1 | genotype), data = bg_nocomp_height)

# Try fitting additional random slopes
height_mod6 <- lmer(height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + elevation_sc +
                    (1+salinity+co2|genotype), data = bg_nocomp_height)

## boundary (singular) fit: see ?isSingular

height_mod7 <- lmer(height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + elevation_sc +
                    (1+salinity+elevation_sc|genotype), data = bg_nocomp_height) 

## boundary (singular) fit: see ?isSingular

height_mod8 <- lmer(height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + elevation_sc +
                    (1+elevation_sc+co2|genotype), data = bg_nocomp_height) # fit
height_mod9 <- lmer(height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + elevation_sc +
                    (1+co2|genotype), data = bg_nocomp_height) # fit
height_mod10 <- lmer(height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + elevation_sc +
                    (1+elevation_sc|genotype), data = bg_nocomp_height) # fit
height_mod11 <- lmer(height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + elevation_sc +
                    (1+salinity|genotype), data = bg_nocomp_height)

## boundary (singular) fit: see ?isSingular

anova(height_mod5, height_mod8)

## refitting model(s) with ML (instead of REML)

## Data: bg_nocomp_height
## Models:
## height_mod5: height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + 
## height_mod5:     elevation_sc + (1 | genotype)
## height_mod8: height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + 
## height_mod8:     elevation_sc + (1 + elevation_sc + co2 | genotype)
##             npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## height_mod5    9 1402.6 1433.4 -692.29   1384.6                     
## height_mod8   14 1406.2 1454.2 -689.12   1378.2 6.3423  5     0.2743

anova(height_mod5, height_mod9) ## This is the best, but not significant

## refitting model(s) with ML (instead of REML)

## Data: bg_nocomp_height
## Models:
## height_mod5: height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + 
## height_mod5:     elevation_sc + (1 | genotype)
## height_mod9: height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + 
## height_mod9:     elevation_sc + (1 + co2 | genotype)
##             npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## height_mod5    9 1402.6 1433.4 -692.29   1384.6                     
## height_mod9   11 1402.0 1439.8 -690.02   1380.0 4.5389  2     0.1034

anova(height_mod5, height_mod10)

## refitting model(s) with ML (instead of REML)

## Data: bg_nocomp_height
## Models:
## height_mod5: height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + 
## height_mod5:     elevation_sc + (1 | genotype)
## height_mod10: height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + 
## height_mod10:     elevation_sc + (1 + elevation_sc | genotype)
##              npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## height_mod5     9 1402.6 1433.4 -692.29   1384.6                    
## height_mod10   11 1405.7 1443.5 -691.87   1383.7 0.836  2     0.6584

# Create plot of fixed effects
Anova(height_mod5) # Just elevation really

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: height_scam_tot
##                     Chisq Df Pr(>Chisq)    
## date_planted_grp   6.9633  1    0.00832 ** 
## age                0.0758  1    0.78304    
## location           0.0986  1    0.75347    
## co2                0.1012  1    0.75038    
## salinity           0.8867  1    0.34638    
## elevation_sc     163.6906  1    < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Refit so it is not scaled for easy plotting
height_modfinal <- lmer(height_scam_tot ~ date_planted_grp + age + location + co2 + salinity + elevation +
                    (1|genotype), data = bg_nocomp_height)

height_means <- emmeans(height_modfinal, ~elevation, at = list(elevation = seq(min(bg_nocomp_height$elevation), max(bg_nocomp_height$elevation), 0.01)), weights = "proportional")

summary(height_means) %>% 
  ggplot(aes(x = elevation, y = emmean)) +
  geom_line(size = 2) +
  geom_line(aes(x = elevation, y = lower.CL), linetype = "dashed")+
  geom_line(aes(x = elevation, y = upper.CL), linetype = "dashed")+
  geom_point(data = bg_nocomp_height, aes(x = elevation, y = height_scam_tot), alpha = 0.2) +
  ylab("Mean stem height (cm)") +
  xlab("Elevation (m NAVD88)") -> k
  
# Random effects
REsim(height_modfinal) %>%
  mutate(groupID = fct_reorder(groupID, mean)) %>% 
  ggplot(aes(x = groupID)) +
  geom_pointrange(aes(ymin = mean - 2*sd, ymax = mean + 2*sd, xmin = groupID, xmax = groupID, y = mean, color = groupID)) +
  coord_flip() +
  ylab("Random effect") +
  xlab("Genotype") +
  geom_hline(aes(yintercept = 0), linetype = "dashed") +
  theme_bw() + scale_x_discrete(position = "top") +
  theme(legend.position = "none") +
  scale_color_manual(values = col_vector) -> l

k + l + plot_layout(widths = c(4,1)) +plot_annotation(tag_levels = 'a') -> height_plot

# Fit null model (no age, location or genotype) and compare to height model
height_mod_null <- lm(height_scam_tot ~ date_planted_grp + co2 + salinity + elevation, data = bg_nocomp_height)
height_null_r2 <- r.squaredGLMM(height_mod_null)[2]
height_r2 <- r.squaredGLMM(height_modfinal)[2]

Average stem width

width_mod <- lmer(width_scam_mid ~ weight_init_sc + date_planted_grp + origin_lab +
                  (age + location + co2 + salinity + elevation_sc)^5 + I(elevation_sc^2) +
                  (1|site_frame) + (1|genotype), data = bg_nocomp)

plot(width_mod)

#step(width_mod, keep = c("co2", "salinity", "elevation_sc", "age", "location"))
summary(width_mod) # site_frame explains nothing

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: width_scam_mid ~ weight_init_sc + date_planted_grp + origin_lab +  
##     (age + location + co2 + salinity + elevation_sc)^5 + I(elevation_sc^2) +  
##     (1 | site_frame) + (1 | genotype)
##    Data: bg_nocomp
## 
## REML criterion at convergence: 362.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6957 -0.5571 -0.0710  0.5351  2.5515 
## 
## Random effects:
##  Groups     Name        Variance  Std.Dev.
##  genotype   (Intercept) 0.0909431 0.30157 
##  site_frame (Intercept) 0.0009338 0.03056 
##  Residual               0.1914315 0.43753 
## Number of obs: 229, groups:  genotype, 31; site_frame, 14
## 
## Fixed effects:
##                                                                       Estimate
## (Intercept)                                                           3.552199
## weight_init_sc                                                        0.067583
## date_planted_grp2                                                    -0.405315
## origin_labutk                                                         0.160306
## agemodern                                                            -0.236577
## locationkirkpatrick                                                  -0.433249
## co2elevated                                                          -0.134332
## salinitysalt                                                         -0.151769
## elevation_sc                                                         -0.469996
## I(elevation_sc^2)                                                     0.006502
## agemodern:locationkirkpatrick                                         0.449555
## agemodern:co2elevated                                                 0.310725
## agemodern:salinitysalt                                                0.338239
## agemodern:elevation_sc                                               -0.072590
## locationkirkpatrick:co2elevated                                       0.206208
## locationkirkpatrick:salinitysalt                                      0.476190
## locationkirkpatrick:elevation_sc                                      0.234863
## co2elevated:salinitysalt                                             -0.042617
## co2elevated:elevation_sc                                              0.045492
## salinitysalt:elevation_sc                                             0.122281
## agemodern:locationkirkpatrick:co2elevated                            -0.122541
## agemodern:locationkirkpatrick:salinitysalt                           -0.369282
## agemodern:locationkirkpatrick:elevation_sc                            0.163139
## agemodern:co2elevated:salinitysalt                                   -0.345161
## agemodern:co2elevated:elevation_sc                                   -0.200902
## agemodern:salinitysalt:elevation_sc                                   0.022314
## locationkirkpatrick:co2elevated:salinitysalt                         -0.067463
## locationkirkpatrick:co2elevated:elevation_sc                         -0.248724
## locationkirkpatrick:salinitysalt:elevation_sc                        -0.189380
## co2elevated:salinitysalt:elevation_sc                                -0.283756
## agemodern:locationkirkpatrick:co2elevated:salinitysalt                0.134748
## agemodern:locationkirkpatrick:co2elevated:elevation_sc               -0.091339
## agemodern:locationkirkpatrick:salinitysalt:elevation_sc              -0.023807
## agemodern:co2elevated:salinitysalt:elevation_sc                       0.445879
## locationkirkpatrick:co2elevated:salinitysalt:elevation_sc             0.616211
## agemodern:locationkirkpatrick:co2elevated:salinitysalt:elevation_sc  -0.262417
##                                                                     Std. Error
## (Intercept)                                                           0.180783
## weight_init_sc                                                        0.034249
## date_planted_grp2                                                     0.199876
## origin_labutk                                                         0.134072
## agemodern                                                             0.194852
## locationkirkpatrick                                                   0.234510
## co2elevated                                                           0.170418
## salinitysalt                                                          0.244546
## elevation_sc                                                          0.107438
## I(elevation_sc^2)                                                     0.039324
## agemodern:locationkirkpatrick                                         0.316275
## agemodern:co2elevated                                                 0.223402
## agemodern:salinitysalt                                                0.194370
## agemodern:elevation_sc                                                0.146213
## locationkirkpatrick:co2elevated                                       0.261749
## locationkirkpatrick:salinitysalt                                      0.255509
## locationkirkpatrick:elevation_sc                                      0.182223
## co2elevated:salinitysalt                                              0.225143
## co2elevated:elevation_sc                                              0.176680
## salinitysalt:elevation_sc                                             0.141817
## agemodern:locationkirkpatrick:co2elevated                             0.350795
## agemodern:locationkirkpatrick:salinitysalt                            0.336443
## agemodern:locationkirkpatrick:elevation_sc                            0.251203
## agemodern:co2elevated:salinitysalt                                    0.304054
## agemodern:co2elevated:elevation_sc                                    0.241305
## agemodern:salinitysalt:elevation_sc                                   0.197682
## locationkirkpatrick:co2elevated:salinitysalt                          0.379913
## locationkirkpatrick:co2elevated:elevation_sc                          0.285121
## locationkirkpatrick:salinitysalt:elevation_sc                         0.266749
## co2elevated:salinitysalt:elevation_sc                                 0.235755
## agemodern:locationkirkpatrick:co2elevated:salinitysalt                0.514829
## agemodern:locationkirkpatrick:co2elevated:elevation_sc                0.408144
## agemodern:locationkirkpatrick:salinitysalt:elevation_sc               0.356388
## agemodern:co2elevated:salinitysalt:elevation_sc                       0.324635
## locationkirkpatrick:co2elevated:salinitysalt:elevation_sc             0.414821
## agemodern:locationkirkpatrick:co2elevated:salinitysalt:elevation_sc   0.591532
##                                                                             df
## (Intercept)                                                          50.726015
## weight_init_sc                                                      157.953703
## date_planted_grp2                                                   179.784008
## origin_labutk                                                        39.580508
## agemodern                                                            54.052246
## locationkirkpatrick                                                  60.502178
## co2elevated                                                          90.361360
## salinitysalt                                                        159.615974
## elevation_sc                                                        168.739610
## I(elevation_sc^2)                                                   171.879127
## agemodern:locationkirkpatrick                                        56.220115
## agemodern:co2elevated                                               174.755364
## agemodern:salinitysalt                                              167.721595
## agemodern:elevation_sc                                              167.284613
## locationkirkpatrick:co2elevated                                     170.293590
## locationkirkpatrick:salinitysalt                                    178.521110
## locationkirkpatrick:elevation_sc                                    174.608314
## co2elevated:salinitysalt                                             76.975497
## co2elevated:elevation_sc                                            168.840090
## salinitysalt:elevation_sc                                           167.422935
## agemodern:locationkirkpatrick:co2elevated                           169.374647
## agemodern:locationkirkpatrick:salinitysalt                          176.998758
## agemodern:locationkirkpatrick:elevation_sc                          179.011320
## agemodern:co2elevated:salinitysalt                                  164.962255
## agemodern:co2elevated:elevation_sc                                  167.803333
## agemodern:salinitysalt:elevation_sc                                 164.677152
## locationkirkpatrick:co2elevated:salinitysalt                        172.679836
## locationkirkpatrick:co2elevated:elevation_sc                        177.411207
## locationkirkpatrick:salinitysalt:elevation_sc                       177.570616
## co2elevated:salinitysalt:elevation_sc                               168.122265
## agemodern:locationkirkpatrick:co2elevated:salinitysalt              169.774274
## agemodern:locationkirkpatrick:co2elevated:elevation_sc              181.543356
## agemodern:locationkirkpatrick:salinitysalt:elevation_sc             179.577296
## agemodern:co2elevated:salinitysalt:elevation_sc                     167.919283
## locationkirkpatrick:co2elevated:salinitysalt:elevation_sc           182.730762
## agemodern:locationkirkpatrick:co2elevated:salinitysalt:elevation_sc 184.191479
##                                                                     t value
## (Intercept)                                                          19.649
## weight_init_sc                                                        1.973
## date_planted_grp2                                                    -2.028
## origin_labutk                                                         1.196
## agemodern                                                            -1.214
## locationkirkpatrick                                                  -1.847
## co2elevated                                                          -0.788
## salinitysalt                                                         -0.621
## elevation_sc                                                         -4.375
## I(elevation_sc^2)                                                     0.165
## agemodern:locationkirkpatrick                                         1.421
## agemodern:co2elevated                                                 1.391
## agemodern:salinitysalt                                                1.740
## agemodern:elevation_sc                                               -0.496
## locationkirkpatrick:co2elevated                                       0.788
## locationkirkpatrick:salinitysalt                                      1.864
## locationkirkpatrick:elevation_sc                                      1.289
## co2elevated:salinitysalt                                             -0.189
## co2elevated:elevation_sc                                              0.257
## salinitysalt:elevation_sc                                             0.862
## agemodern:locationkirkpatrick:co2elevated                            -0.349
## agemodern:locationkirkpatrick:salinitysalt                           -1.098
## agemodern:locationkirkpatrick:elevation_sc                            0.649
## agemodern:co2elevated:salinitysalt                                   -1.135
## agemodern:co2elevated:elevation_sc                                   -0.833
## agemodern:salinitysalt:elevation_sc                                   0.113
## locationkirkpatrick:co2elevated:salinitysalt                         -0.178
## locationkirkpatrick:co2elevated:elevation_sc                         -0.872
## locationkirkpatrick:salinitysalt:elevation_sc                        -0.710
## co2elevated:salinitysalt:elevation_sc                                -1.204
## agemodern:locationkirkpatrick:co2elevated:salinitysalt                0.262
## agemodern:locationkirkpatrick:co2elevated:elevation_sc               -0.224
## agemodern:locationkirkpatrick:salinitysalt:elevation_sc              -0.067
## agemodern:co2elevated:salinitysalt:elevation_sc                       1.373
## locationkirkpatrick:co2elevated:salinitysalt:elevation_sc             1.485
## agemodern:locationkirkpatrick:co2elevated:salinitysalt:elevation_sc  -0.444
##                                                                     Pr(>|t|)
## (Intercept)                                                          < 2e-16
## weight_init_sc                                                        0.0502
## date_planted_grp2                                                     0.0441
## origin_labutk                                                         0.2389
## agemodern                                                             0.2300
## locationkirkpatrick                                                   0.0696
## co2elevated                                                           0.4326
## salinitysalt                                                          0.5357
## elevation_sc                                                        2.12e-05
## I(elevation_sc^2)                                                     0.8689
## agemodern:locationkirkpatrick                                         0.1607
## agemodern:co2elevated                                                 0.1660
## agemodern:salinitysalt                                                0.0837
## agemodern:elevation_sc                                                0.6202
## locationkirkpatrick:co2elevated                                       0.4319
## locationkirkpatrick:salinitysalt                                      0.0640
## locationkirkpatrick:elevation_sc                                      0.1991
## co2elevated:salinitysalt                                              0.8504
## co2elevated:elevation_sc                                              0.7971
## salinitysalt:elevation_sc                                             0.3898
## agemodern:locationkirkpatrick:co2elevated                             0.7273
## agemodern:locationkirkpatrick:salinitysalt                            0.2739
## agemodern:locationkirkpatrick:elevation_sc                            0.5169
## agemodern:co2elevated:salinitysalt                                    0.2579
## agemodern:co2elevated:elevation_sc                                    0.4063
## agemodern:salinitysalt:elevation_sc                                   0.9103
## locationkirkpatrick:co2elevated:salinitysalt                          0.8593
## locationkirkpatrick:co2elevated:elevation_sc                          0.3842
## locationkirkpatrick:salinitysalt:elevation_sc                         0.4787
## co2elevated:salinitysalt:elevation_sc                                 0.2304
## agemodern:locationkirkpatrick:co2elevated:salinitysalt                0.7938
## agemodern:locationkirkpatrick:co2elevated:elevation_sc                0.8232
## agemodern:locationkirkpatrick:salinitysalt:elevation_sc               0.9468
## agemodern:co2elevated:salinitysalt:elevation_sc                       0.1714
## locationkirkpatrick:co2elevated:salinitysalt:elevation_sc             0.1391
## agemodern:locationkirkpatrick:co2elevated:salinitysalt:elevation_sc   0.6578
##                                                                        
## (Intercept)                                                         ***
## weight_init_sc                                                      .  
## date_planted_grp2                                                   *  
## origin_labutk                                                          
## agemodern                                                              
## locationkirkpatrick                                                 .  
## co2elevated                                                            
## salinitysalt                                                           
## elevation_sc                                                        ***
## I(elevation_sc^2)                                                      
## agemodern:locationkirkpatrick                                          
## agemodern:co2elevated                                                  
## agemodern:salinitysalt                                              .  
## agemodern:elevation_sc                                                 
## locationkirkpatrick:co2elevated                                        
## locationkirkpatrick:salinitysalt                                    .  
## locationkirkpatrick:elevation_sc                                       
## co2elevated:salinitysalt                                               
## co2elevated:elevation_sc                                               
## salinitysalt:elevation_sc                                              
## agemodern:locationkirkpatrick:co2elevated                              
## agemodern:locationkirkpatrick:salinitysalt                             
## agemodern:locationkirkpatrick:elevation_sc                             
## agemodern:co2elevated:salinitysalt                                     
## agemodern:co2elevated:elevation_sc                                     
## agemodern:salinitysalt:elevation_sc                                    
## locationkirkpatrick:co2elevated:salinitysalt                           
## locationkirkpatrick:co2elevated:elevation_sc                           
## locationkirkpatrick:salinitysalt:elevation_sc                          
## co2elevated:salinitysalt:elevation_sc                                  
## agemodern:locationkirkpatrick:co2elevated:salinitysalt                 
## agemodern:locationkirkpatrick:co2elevated:elevation_sc                 
## agemodern:locationkirkpatrick:salinitysalt:elevation_sc                
## agemodern:co2elevated:salinitysalt:elevation_sc                        
## locationkirkpatrick:co2elevated:salinitysalt:elevation_sc              
## agemodern:locationkirkpatrick:co2elevated:salinitysalt:elevation_sc    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 36 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

width_mod2 <- lmer(width_scam_mid ~ weight_init_sc + date_planted_grp + origin_lab +
                  (age + location + co2 + salinity + elevation_sc)^5 + I(elevation_sc^2) +
                  (1|genotype), data = bg_nocomp)

step(width_mod2, keep = c("co2", "salinity", "elevation_sc", "age", "location"))

## Backward reduced random-effect table:
## 
##                Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                      38 -181.14 438.28                         
## (1 | genotype)          0   37 -198.30 470.60 34.323  1  4.668e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                                        Eliminated  Sum Sq Mean Sq NumDF   DenDF
## I(elevation_sc^2)                               1 0.00558 0.00558     1 176.642
## age:location:co2:salinity:elevation_sc          2 0.03915 0.03915     1 190.963
## age:location:co2:salinity                       3 0.00722 0.00722     1 175.156
## age:location:salinity:elevation_sc              4 0.04394 0.04394     1 178.924
## age:location:co2:elevation_sc                   5 0.12646 0.12646     1 184.086
## age:location:co2                                6 0.01269 0.01269     1 182.208
## age:location:elevation_sc                       7 0.03928 0.03928     1 180.835
## age:location:salinity                           8 0.28527 0.28527     1 185.334
## age:location                                    9 0.21127 0.21127     1  26.253
## origin_lab                                     10 0.35404 0.35404     1  43.470
## age:co2:salinity:elevation_sc                  11 0.38770 0.38770     1 191.499
## age:co2:elevation_sc                           12 0.02575 0.02575     1 187.509
## age:salinity:elevation_sc                      13 0.22189 0.22189     1 183.576
## age:elevation_sc                               14 0.02047 0.02047     1 185.188
## age:co2:salinity                               15 0.25857 0.25857     1 183.807
## age:co2                                        16 0.12120 0.12120     1 191.546
## age:salinity                                   17 0.16924 0.16924     1 189.549
## location:co2:salinity:elevation_sc             18 0.39122 0.39122     1 205.888
## location:salinity:elevation_sc                 19 0.00054 0.00054     1 190.663
## location:co2:elevation_sc                      20 0.00836 0.00836     1 193.646
## location:co2:salinity                          21 0.01331 0.01331     1 190.110
## co2:salinity:elevation_sc                      22 0.08143 0.08143     1 201.648
## location:co2                                   23 0.14072 0.14072     1 196.758
## co2:elevation_sc                               24 0.40056 0.40056     1 197.858
## weight_init_sc                                 25 0.61674 0.61674     1 204.488
## salinity:elevation_sc                          26 0.60474 0.60474     1 197.007
## co2:salinity                                   27 0.66749 0.66749     1 195.130
## date_planted_grp                                0 1.61835 1.61835     1 201.850
## age                                             0 0.06275 0.06275     1  27.501
## co2                                             0 0.00276 0.00276     1 203.339
## location:salinity                               0 0.81352 0.81352     1 199.950
## location:elevation_sc                           0 1.48509 1.48509     1 197.706
##                                        F value   Pr(>F)   
## I(elevation_sc^2)                       0.0290 0.864957   
## age:location:co2:salinity:elevation_sc  0.2048 0.651390   
## age:location:co2:salinity               0.0379 0.845873   
## age:location:salinity:elevation_sc      0.2317 0.630840   
## age:location:co2:elevation_sc           0.6701 0.414084   
## age:location:co2                        0.0675 0.795378   
## age:location:elevation_sc               0.2099 0.647398   
## age:location:salinity                   1.5321 0.217359   
## age:location                            1.1276 0.297966   
## origin_lab                              1.8885 0.176419   
## age:co2:salinity:elevation_sc           2.0565 0.153194   
## age:co2:elevation_sc                    0.1361 0.712575   
## age:salinity:elevation_sc               1.1785 0.279079   
## age:elevation_sc                        0.1086 0.742070   
## age:co2:salinity                        1.3792 0.241763   
## age:co2                                 0.6466 0.422312   
## age:salinity                            0.9042 0.342882   
## location:co2:salinity:elevation_sc      2.0909 0.149695   
## location:salinity:elevation_sc          0.0029 0.957164   
## location:co2:elevation_sc               0.0447 0.832776   
## location:co2:salinity                   0.0716 0.789320   
## co2:salinity:elevation_sc               0.4402 0.507804   
## location:co2                            0.7630 0.383471   
## co2:elevation_sc                        2.1761 0.141757   
## weight_init_sc                          3.3336 0.069335 . 
## salinity:elevation_sc                   3.2392 0.073427 . 
## co2:salinity                            3.5374 0.061490 . 
## date_planted_grp                        8.4408 0.004078 **
## age                                     0.3273 0.571898   
## co2                                     0.0144 0.904573   
## location:salinity                       4.2430 0.040707 * 
## location:elevation_sc                   7.7457 0.005906 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## width_scam_mid ~ date_planted_grp + age + location + co2 + salinity + 
##     elevation_sc + (1 | genotype) + location:salinity + location:elevation_sc

# Fit proposed model
width_mod3 <- lmer(width_scam_mid ~ date_planted_grp + age + co2 +
                     location*salinity + location*elevation_sc + (1|genotype), data = bg_nocomp)

# Diagnostic plots
plot_model(width_mod3, type = "diag") # Looks good

## [[1]]

## `geom_smooth()` using formula 'y ~ x'

## 
## [[2]]
## [[2]]$genotype

## `geom_smooth()` using formula 'y ~ x'

## 
## 
## [[3]]

## 
## [[4]]

## `geom_smooth()` using formula 'y ~ x'

# See if we can add any random effects (all additive)

width_mod4 <- lmer(width_scam_mid ~ date_planted_grp + age + co2 +
                     location*salinity + location*elevation_sc +
                     (1+co2+salinity|genotype), data = bg_nocomp) # fits
width_mod5 <- lmer(width_scam_mid ~ date_planted_grp + age + co2 +
                     location*salinity + location*elevation_sc +
                     (1+co2+elevation_sc|genotype), data = bg_nocomp)

## boundary (singular) fit: see ?isSingular

width_mod6 <- lmer(width_scam_mid ~ date_planted_grp + age + co2 +
                     location*salinity + location*elevation_sc +
                     (1+salinity+elevation_sc|genotype), data = bg_nocomp) # fits
width_mod7 <- lmer(width_scam_mid ~ date_planted_grp + age + co2 +
                     location*salinity + location*elevation_sc +
                     (1+salinity|genotype), data = bg_nocomp) # fits
width_mod8 <- lmer(width_scam_mid ~ date_planted_grp + age + co2 +
                     location*salinity + location*elevation_sc +
                     (1+elevation_sc|genotype), data = bg_nocomp) # fits
width_mod9 <- lmer(width_scam_mid ~ date_planted_grp + age + co2 +
                     location*salinity + location*elevation_sc +
                     (1+co2|genotype), data = bg_nocomp)

## boundary (singular) fit: see ?isSingular

anova(width_mod3, width_mod4)

## refitting model(s) with ML (instead of REML)

## Data: bg_nocomp
## Models:
## width_mod3: width_scam_mid ~ date_planted_grp + age + co2 + location * salinity + 
## width_mod3:     location * elevation_sc + (1 | genotype)
## width_mod4: width_scam_mid ~ date_planted_grp + age + co2 + location * salinity + 
## width_mod4:     location * elevation_sc + (1 + co2 + salinity | genotype)
##            npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## width_mod3   11 331.06 368.83 -154.53   309.06                     
## width_mod4   16 337.85 392.79 -152.93   305.86 3.2055  5     0.6683

anova(width_mod3, width_mod6)

## refitting model(s) with ML (instead of REML)

## Data: bg_nocomp
## Models:
## width_mod3: width_scam_mid ~ date_planted_grp + age + co2 + location * salinity + 
## width_mod3:     location * elevation_sc + (1 | genotype)
## width_mod6: width_scam_mid ~ date_planted_grp + age + co2 + location * salinity + 
## width_mod6:     location * elevation_sc + (1 + salinity + elevation_sc | 
## width_mod6:     genotype)
##            npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## width_mod3   11 331.06 368.83 -154.53   309.06                     
## width_mod6   16 336.74 391.68 -152.37   304.74 4.3164  5     0.5048

anova(width_mod3, width_mod7)

## refitting model(s) with ML (instead of REML)

## Data: bg_nocomp
## Models:
## width_mod3: width_scam_mid ~ date_planted_grp + age + co2 + location * salinity + 
## width_mod3:     location * elevation_sc + (1 | genotype)
## width_mod7: width_scam_mid ~ date_planted_grp + age + co2 + location * salinity + 
## width_mod7:     location * elevation_sc + (1 + salinity | genotype)
##            npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## width_mod3   11 331.06 368.83 -154.53   309.06                     
## width_mod7   13 332.58 377.22 -153.29   306.58 2.4836  2     0.2889

anova(width_mod3, width_mod8)

## refitting model(s) with ML (instead of REML)

## Data: bg_nocomp
## Models:
## width_mod3: width_scam_mid ~ date_planted_grp + age + co2 + location * salinity + 
## width_mod3:     location * elevation_sc + (1 | genotype)
## width_mod8: width_scam_mid ~ date_planted_grp + age + co2 + location * salinity + 
## width_mod8:     location * elevation_sc + (1 + elevation_sc | genotype)
##            npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## width_mod3   11 331.06 368.83 -154.53   309.06                     
## width_mod8   13 333.72 378.36 -153.86   307.72 1.3375  2     0.5123

# None are significant

# Plot the fixed effects
Anova(width_mod3)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: width_scam_mid
##                          Chisq Df Pr(>Chisq)    
## date_planted_grp        8.4408  1   0.003669 ** 
## age                     0.3273  1   0.567251    
## co2                     0.0144  1   0.904454    
## location                0.0168  1   0.896926    
## salinity                1.1383  1   0.286019    
## elevation_sc          181.7487  1  < 2.2e-16 ***
## location:salinity       4.2430  1   0.039412 *  
## location:elevation_sc   7.7457  1   0.005384 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Make fixed effect plot for location by elevation
width_means <- emmeans(width_mod3, ~location:elevation_sc,
                    at= list(elevation_sc = seq(-1.72, 1.54, 0.1)), type = "response")

width_means_out <- summary(width_means)

bg_nocomp_plot %>% 
  mutate(location = case_when(location == "corn" ~ "Corn Island",
                              T ~ "Kirkpatrick Marsh")) -> bg_nocomp_plot

width_means_out %>% 
  mutate(location = case_when(location == "corn" ~ "Corn Island",
                              T ~ "Kirkpatrick Marsh")) %>% 
  mutate(elevation = elevation_sc*scaling + centering) %>% 
  ggplot(aes(x = elevation, y = emmean, color = location)) +
  geom_line(size = 1.5) +
  geom_point(data = bg_nocomp_plot, aes(x = elevation, y = width_scam_mid), alpha = 0.5) +
  theme_bw() + scale_color_manual(values = c("#1b9e77", "#7570b3")) +
  ylab("Mean stem width (mm)") +
  xlab("Elevation (m NAVD88)") +
  guides(color=guide_legend("")) +
  theme(legend.position = "top") -> m
  
# Same but for location:salinity interaction
width_mod3 <- lmer(width_scam_mid ~ date_planted_grp + age + co2 +
                     location*salinity + location*elevation_sc + (1|genotype), data = bg_nocomp_plot)

emmip(width_mod3, location ~ salinity, CI = T, weights = "proportional") +
  geom_jitter(data = bg_nocomp_rs_plot, aes(x = salinity, y = width_scam_mid, color = location),
              height = 0, width = 0.1, alpha = 0.3, shape = 1) +
  ylab("Mean stem width (mm)") +
  scale_color_manual(values = c("#1b9e77", "#7570b3")) + theme_bw() +
  theme(legend.position = "top") + guides(color=guide_legend("")) +
  scale_size_manual(values=c(2,2)) -> n
 

# Random effects plot
REsim(width_mod3) %>%
  mutate(groupID = fct_reorder(groupID, mean)) %>% 
  ggplot(aes(x = groupID)) +
  geom_pointrange(aes(ymin = mean - 2*sd, ymax = mean + 2*sd, xmin = groupID, xmax = groupID, y = mean, color = groupID)) +
  coord_flip() +
  ylab("Random effect") +
  xlab("Genotype") +
  geom_hline(aes(yintercept = 0), linetype = "dashed") +
  theme_bw() + scale_x_discrete(position = "top") +
  theme(legend.position = "none") +
  scale_color_manual(values = col_vector) -> o

mn <- cowplot::plot_grid(m, n, ncol = 1, labels = c("a", "b"))

## Warning: Removed 1 rows containing missing values (geom_point).

o_alone <- cowplot::plot_grid(o, labels = c("c"))
mn+o_alone+plot_layout(widths = c(4,1)) -> width_plot

# Fit a model without age, location or genotype
width_mod_null <- lm(width_scam_mid ~ date_planted_grp + co2 +
                     salinity + elevation_sc, data = bg_nocomp_plot)

width_null_r2 <- r.squaredGLMM(width_mod_null)[2]
width_r2 <- r.squaredGLMM(width_mod3)[2]

Belowground biomass

bgb_mod <- lmer(total_bg ~ weight_init_sc + date_planted_grp + origin_lab +
                  (age + location + co2 + salinity + elevation_sc)^5 + I(elevation_sc^2) +
                  (1|site_frame) + (1|genotype), data = bg_nocomp)

#step(bgb_mod)
VarCorr(bgb_mod)

##  Groups     Name        Std.Dev.
##  genotype   (Intercept) 1.40014 
##  site_frame (Intercept) 0.27376 
##  Residual               3.38998

bgb_mod2 <- lmer(total_bg ~ weight_init_sc + date_planted_grp + origin_lab +
                  (age + location + co2 + salinity + elevation_sc)^5 + I(elevation_sc^2) +
                  (1|genotype), data = bg_nocomp)

step(bgb_mod2, keep = c("co2", "elevation_sc", "salinity", "age", "location"))

## Backward reduced random-effect table:
## 
##                Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)   
## <none>                      38 -565.32 1206.7                        
## (1 | genotype)          0   37 -568.81 1211.6 6.9719  1    0.00828 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                                        Eliminated Sum Sq Mean Sq NumDF   DenDF
## age:location:co2:salinity:elevation_sc          1   0.30    0.30     1 191.632
## origin_lab                                      2   0.80    0.80     1  31.107
## age:location:co2:elevation_sc                   3   2.13    2.13     1 183.313
## age:location:co2:salinity                       4   6.42    6.42     1 176.485
## age:location:salinity:elevation_sc              5   8.66    8.66     1 180.800
## age:location:elevation_sc                       6   0.00    0.00     1 181.213
## age:location:salinity                           7   3.31    3.31     1 185.205
## age:co2:salinity:elevation_sc                   8  13.57   13.57     1 196.165
## age:co2:elevation_sc                            9  13.05   13.05     1 190.620
## location:co2:salinity:elevation_sc             10  13.26   13.26     1 200.299
## co2:salinity:elevation_sc                      11   0.65    0.65     1 200.325
## location:co2:elevation_sc                      12   2.49    2.49     1 194.922
## location:co2:salinity                          13   2.76    2.76     1 185.776
## age:salinity:elevation_sc                      14  11.66   11.66     1 190.162
## age:elevation_sc                               15   6.95    6.95     1 190.408
## location:salinity:elevation_sc                 16  11.15   11.15     1 193.654
## location:elevation_sc                          17   0.10    0.10     1 191.214
## location:salinity                              18   2.30    2.30     1 196.923
## co2:elevation_sc                               19  12.25   12.25     1 200.383
## salinity:elevation_sc                          20  24.80   24.80     1 197.169
## age:location:co2                               21  24.55   24.55     1 202.999
## location:co2                                   22   1.64    1.64     1 204.121
## age:location                                   23  11.33   11.33     1  25.241
## weight_init_sc                                  0  53.05   53.05     1 213.800
## date_planted_grp                                0  53.89   53.89     1 205.811
## location                                        0   1.50    1.50     1  27.270
## elevation_sc                                    0   0.02    0.02     1 199.343
## I(elevation_sc^2)                               0 386.08  386.08     1 204.957
## age:co2:salinity                                0  45.46   45.46     1 195.963
##                                        F value    Pr(>F)    
## age:location:co2:salinity:elevation_sc  0.0260   0.87200    
## origin_lab                              0.0699   0.79323    
## age:location:co2:elevation_sc           0.1854   0.66726    
## age:location:co2:salinity               0.5606   0.45500    
## age:location:salinity:elevation_sc      0.7585   0.38496    
## age:location:elevation_sc               0.0000   0.99442    
## age:location:salinity                   0.2915   0.58994    
## age:co2:salinity:elevation_sc           1.1995   0.27476    
## age:co2:elevation_sc                    1.1492   0.28508    
## location:co2:salinity:elevation_sc      1.1603   0.28270    
## co2:salinity:elevation_sc               0.0563   0.81261    
## location:co2:elevation_sc               0.2169   0.64193    
## location:co2:salinity                   0.2425   0.62296    
## age:salinity:elevation_sc               1.0274   0.31206    
## age:elevation_sc                        0.6120   0.43500    
## location:salinity:elevation_sc          0.9874   0.32161    
## location:elevation_sc                   0.0093   0.92346    
## location:salinity                       0.2048   0.65138    
## co2:elevation_sc                        1.0961   0.29637    
## salinity:elevation_sc                   2.2101   0.13871    
## age:location:co2                        2.1710   0.14218    
## location:co2                            0.1442   0.70451    
## age:location                            0.9984   0.32718    
## weight_init_sc                          4.6787   0.03165 *  
## date_planted_grp                        4.7525   0.03039 *  
## location                                0.1321   0.71908    
## elevation_sc                            0.0021   0.96312    
## I(elevation_sc^2)                      34.0485 2.075e-08 ***
## age:co2:salinity                        4.0094   0.04663 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## total_bg ~ weight_init_sc + date_planted_grp + age + location + 
##     co2 + salinity + elevation_sc + I(elevation_sc^2) + (1 | 
##     genotype) + age:co2 + age:salinity + co2:salinity + age:co2:salinity

bgb_mod3 <- lmer(total_bg ~ weight_init_sc + date_planted_grp + age + co2 + salinity + elevation +
    + location + I(elevation^2) + (1 | genotype) + age*co2*salinity, data = bg_nocomp_plot)

# Diagnostics
plot_model(bgb_mod3, type = "diag") # Pretty!

## [[1]]

## `geom_smooth()` using formula 'y ~ x'

## 
## [[2]]
## [[2]]$genotype

## `geom_smooth()` using formula 'y ~ x'

## 
## 
## [[3]]

## 
## [[4]]

## `geom_smooth()` using formula 'y ~ x'

# Fixed effects plots
car::Anova(bgb_mod3)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: total_bg
##                    Chisq Df Pr(>Chisq)    
## weight_init_sc    4.6787  1   0.030539 *  
## date_planted_grp  4.7525  1   0.029256 *  
## age               0.0091  1   0.923897    
## co2              13.8404  1   0.000199 ***
## salinity          2.6561  1   0.103153    
## elevation        32.5434  1  1.166e-08 ***
## location          0.1321  1   0.716282    
## I(elevation^2)   34.0485  1  5.376e-09 ***
## age:co2           0.0088  1   0.925363    
## age:salinity      0.2140  1   0.643665    
## co2:salinity      2.8146  1   0.093409 .  
## age:co2:salinity  4.0094  1   0.045247 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

bgb_means <- emmeans::emmeans(bgb_mod3, ~co2:salinity:age:elevation, at= list(elevation = seq(0.15, 0.55, 0.01)), type = "response")

bgb_means_out <- summary(bgb_means)

bgb_means_out %>% 
  ggplot(aes(x = elevation, y = emmean, color = co2)) + 
  geom_line(size = 1.5, alpha= 0.8) +
  geom_line(aes(x = elevation, y = lower.CL), linetype = "dashed") + 
  geom_line(aes(x = elevation, y = upper.CL), linetype = "dashed") + 
  facet_grid(age~salinity) +
  geom_point(data = bg_nocomp_plot, aes(x = elevation, y = total_bg), alpha = 0.5) +
  theme_bw() +
  scale_color_manual(values = c("gray47", "#66a61e")) +
  ylab("Belowground biomass (g)") +
  xlab("Elevation (m NAVD88)") +
  guides(color=guide_legend(title=expression(CO[2]))) +
  theme(legend.position = "top") -> p

# See if we can add random effects
bgb_mod4 <- lmer(total_bg ~ weight_init_sc + date_planted_grp + age + co2 + salinity + elevation_sc + location + I(elevation_sc^2) + (1 + salinity + elevation_sc | genotype) + age*co2*salinity, data = bg_nocomp_plot)

## boundary (singular) fit: see ?isSingular

bgb_mod5 <- lmer(total_bg ~ weight_init_sc + date_planted_grp + age + co2 + salinity + elevation_sc + location + I(elevation_sc^2) + (1 + co2 + elevation_sc | genotype) + age*co2*salinity, data = bg_nocomp_plot)

## boundary (singular) fit: see ?isSingular

bgb_mod6 <- lmer(total_bg ~ weight_init_sc + date_planted_grp + age + co2 + salinity + elevation_sc + location + I(elevation_sc^2) + (1 + salinity + co2 | genotype) + age*co2*salinity, data = bg_nocomp_plot)

## boundary (singular) fit: see ?isSingular

bgb_mod7 <- lmer(total_bg ~ weight_init_sc + date_planted_grp + age + co2 + salinity + elevation_sc + location + I(elevation_sc^2) + (1 + co2 | genotype) + age*co2*salinity, data = bg_nocomp_plot)

## boundary (singular) fit: see ?isSingular

bgb_mod8 <- lmer(total_bg ~ weight_init_sc + date_planted_grp + age + co2 + salinity + elevation_sc + location + I(elevation_sc^2) + (1 + salinity | genotype) + age*co2*salinity, data = bg_nocomp_plot)

## boundary (singular) fit: see ?isSingular

bgb_mod9 <- lmer(total_bg ~ weight_init_sc + date_planted_grp + age + co2 + salinity + elevation_sc + location + I(elevation_sc^2) + (1 + elevation_sc | genotype) + age*co2*salinity, data = bg_nocomp_plot)

## boundary (singular) fit: see ?isSingular

# No random slopes

# Create caterpillar graph
REsim(bgb_mod3) %>%
  mutate(groupID = fct_reorder(groupID, mean)) %>% 
  ggplot(aes(x = groupID)) +
  geom_pointrange(aes(ymin = mean - 2*sd, ymax = mean + 2*sd, xmin = groupID, xmax = groupID, y = mean, color = groupID)) +
  coord_flip() +
  ylab("Random effect") +
  xlab("Genotype") +
  geom_hline(aes(yintercept = 0), linetype = "dashed") +
  theme_bw() + scale_x_discrete(position = "top") +
  theme(legend.position = "none") +
  scale_color_manual(values = col_vector)-> q

p + q + plot_layout(widths = c(4,1)) -> bgb_plot

# Build null model
bgb_null <- lm(total_bg ~ weight_init_sc + date_planted_grp + elevation +
     + I(elevation^2) + co2*salinity, data = bg_nocomp_plot)
bg_null_r2 <- r.squaredGLMM(bgb_null)[2]
bg_r2 <- r.squaredGLMM(bgb_mod3)[2]

## Warning: Removed 2 rows containing missing values (geom_point).

## Warning: Removed 3 rows containing missing values (geom_point).

## Warning: Removed 3 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

trait	null r2	evo r2
aboveground biomass	0.35	0.46
root-to-shoot ratio	0.60	0.69
belowground biomass	0.23	0.32
root distribution parameter	0.46	0.76
mean stem height	0.38	0.60
mean stem width	0.42	0.62